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PREDICTION  OF  TRANSITORY  STALL  IN  TWO-DIMENSIONAL  DIFFUSERS 

A method  has  been  developed  that  predicts  the  performance  of  diffu- 
sers operating  in  the  transitory  stall  mode  of  the  flow  regime  chart  of 
Fox  and  Kline.  The  calculations  are  accurate  within  ± 6%,  which  is  of 
the  same  order  as  the  uncertainty  in  the  data  for  diffusers  with  diver- 
gence angles  that  are  1.2  times  that  at  which  line  a-a  occurs.  This 
corresponds  approximately  to  the  line  of  appreciable  stall. 

Singular  behavior  in  the  neighborhood  of  detachment  is  avoided  by 
simultaneous  calculation  of  the  inviscid  core  and  the  shear  layers.  A 
new  boundary  layer  scheme  using  Bradshaw's  entrainment-maximum  shear  cor- 
relation is  developed  that  is  valid  for  both  attached  and  detached  flows. 

The  irrotational  core  is  first  assumed  one-dimensional  and  then  extended 
to  the  two-dimensional  case  by  an  iterative  scheme  consisting  of  alter- 
nate calculations  of  the  boundary  layers  and  Laplace's  equation  in  the 
core . 

The  basic  boundary  layer  method  is  shown  to  be  of  comparable  accu- 
racy as  the  best  calculation  presented  in  the  1968  Stanford  Conference 
on  Computation  of  Turbulent  Boundary  Layers.  When  compared  against  the 
data  maps  of  Reneau  et  al.  and  the  measurements  of  Carlson  et  al.,  the 
one-dimensional  core  model  gives  excellent  agreement  for  the  streamwise 
distribution  of  the  shape  factor  H,  the  displacement  thickness,  and 
skin  friction  coefficient  C^/2,  as  well  as  for  the  locations  of  inter- 
mittent detachment  and  time-averaged  zero  wall  shear.  The  two-dimensional 
model  predicts  the  same  quantities  to  the  accuracy  in  the  data  for  the 
flow  of  Strickland  and  Simpson.  However,  in  the  reversed  flow  portion, 
the  predicted  skin  friction  is  somewhat  low,  and  the  entrainment  much 
too  high.  In  all  cases,  the  largest  deviation  from  data  occurs  in  the 
region  between  intermittent  detachment  and  the  location  of  time-averaged 
zero  wall  shear.  Complete  verification  of  the  method,  or  its  improvement 
in  this  zone,  must  await  further  data. 
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CHAPTER  ONE 
INTRODUCTION 


A.  Ob j ective 

The  objective  of  this  investigation  is  to  develop  a method  for  pre- 
dicting the  performance  of  two-dimensional  (2-D)  diffusers  operating  in 
the  "unstalled"  and  "transitory  stalled"  regimes  of  the  diffuser  perfor- 
mance chart  of  Fox  and  Kline  [1],  as  shown  in  Fig.  1. 

A typical  curve  of  static  pressure  recovery,  Cp,  as  a function  of 
the  divergence  angle  20  is  shown  in  Fig.  A.^  Line  a-a  on  both  fig- 
ures represents  the  approximate  dividing  line  between  the  unstalled  and 
the  transitory  stalled  regimes.  Current  calculation  methods  [4,5]  can 
make  successful  predictions  in  the  shaded  zones  only.  These  consist  of 
the  fully  stalled  regime  and  the  unstalled  zone  for  diffusers  with 
20/20  = 0.6  to  0.8.  It  will  be  noted  that  the  peak  pressure  recovery, 

•k 

Cp  , occurs  in  the  transitory  stall  regime,  where  the  flow  is  unsteady 
and  the  boundary  layers  (b.l.'s)  along  the  diffuser  walls  are  partially 
detached.  (Separated  or  stalled  b.l.'s  shall  be  referred  to  as  detached 
b.l.'s  to  avoid  confusion  between  stalled  b.l.'s  and  stalled  diffusers.) 

k 

The  ability  to  predict  diffuser  performance  in  the  region  near  Cp 
is  of  obvious  interest  to  the  designer  of  flow  equipment.  However,  a 
prerequisite  to  being  able  to  do  this  is  the  calculation  of  b.l.'s  which 
are  attached  and  partially  detached.  Accordingly,  the  first  few  chapters 
of  this  report  are  devoted  to  the  development  of  such  a turbulent  bound- 
ary layer  prediction  method  (TBLPM) , which  is  then  used  to  predict  diffu- 

* 

sers  operating  in  the  region  near  Cp  . 

B.  Cyclic  Iteration 

The  classical  method  for  predicting  the  development  of  b.l.'s  is  to 
prescribe  the  pressure  gradient  dp/dx  in  the  flow  direction  and  calcu- 
late the  dependent  b.l.  parameters  through  a parabolic  marching  scheme 

Figures  with  lettered  titles  (Fig.  A,  Fig.  B,  etc.)  are  embodied 
within  the  text.  Numbered  figures  (Fig.  1,  Fig.  2,  etc.)  are  collected 
at  the  end. 
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Fig.  A.  Behavior  of  C with  increasing  area  ratio  or  2 A. 


[2,3].  The  a priori  imposition  of  the  pressure  gradient  implies  that  the 
b.l.  does  not  greatly  affect  the  free  stream  velocity,  an  assumption  that 
is  true  only  for  attached  flows  with  b.l.'s  that  are  thin  compared  to 
passage  height. 

The  pressure  gradient  acting  on  the  b.l.  is  in  fact  the  result  of 
mutual  interaction  between  itself  and  the  adjacent  irrotational  fluid. 

This  interaction  assumes  an  increasingly  important  role  in  adverse  pres- 
sure gradients,  as  the  "blockage"  of  the  b.l.  becomes  greater  and  greater. 
Finally,  for  flows  at  and  near  detachment,  it  will  be  shown  to  be  the 
controlling  factor  in  determining  whether  or  not  such  calculations  can 
be  made  to  converge. 

This  problem  is  much  more  serious  for  internal  flow  than  for  exter- 
nal flow,  because  in  internal  flow  the  irrotational  core  is  confined  be- 
tween b.l.'s  growing  on  the  bounding  walls  and  the  blockage  is  large 
enough  in  typical  passages  to  cause  a substantial  amount  of  mutual  inter- 
action. 

Several  schemes  for  fully  stalled  and  attached  flows  have  recently 
appeared  [4,5],  wherein  the  classical  turbulent  boundary  layer  prediction 
methods  (TBLPM's)  with  prescribed  pressure  gradient  have  been  coupled 
with  an  inviscid  core  and  the  calculation  iterated  to  closure  using  a 
scheme  such  as  shown  in  Fig.  B.  An  initial  estimate  of  the  pressure  gra- 
dient is  impressed  upon  and  used  to  calculate  the  b.l.'s,  which  in  turn 

* 

supply  an  estimate  of  the  displacement  thickness,  6 . The  blockage  is 
subtracted  from  the  channel  width,  giving  a new  body  shape,  which  is  then 
used  to  calculate  the  new  pressure  gradient  and  so  on,  hopefully  to  con- 
vergence. We  shall  call  such  schemes  "CYCLIC  ITERATION". 

Consider  the  case  of  a rapidly  detaching  flow,  such  as  in  a stalled 

* 

diffuser.  As  the  b.l.  approaches  detachment,  6 grows  very  rapidly. 

In  a real  physical  situation,  this  rapid  increase  in  blockage  will  result 
in  a simultaneous  decrease  in  pressure  gradient.  This  is  so  because  the 
mutual  interaction  between  the  b.l.  and  the  potential  core  decreases  the 
effective  flow  channel  (EFC)  available  to  the  outer  flow,  thereby  relax- 
ing the  pressure  gradient  as  shown  in  Fig.  C. 

In  cyclic  iteration,  however,  the  pressure  gradient  is  fixed  before- 
hand for  the  entire  iteration,  and  there  is  no  mechanism  available  to 


3 


DIFFUSER 
GEOMETRY 
(ZERO  BLOCKAGE) 


4 


Fig.  B.  Flowchart  of  cyclic  iteration  as  applied  to  the  calculation 
of  diffusers. 
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Fig.  C.  Behavior  of  6 in  a real  flow. 


* 

reduce  the  pressure  gradient  in  reaction  to  the  sudden  growth  of  6 

(Fig.  D)  in  that  iteration.  The  adverse  pressure  gradient  is  therefore 

•k 

maintained  unchanged,  resulting  in  runaway  growth  of  6 and  catas- 
trophic failure  of  the  prediction  scheme.  This  is  the  so-called  "sepa- 
ration singularity",  the  effects  of  which  can  be  seen  rather  dramatically 
in  many  of  the  "separating  flows"  of  reference  [3]. 

A related  effect  is  the  inability  of  the  calculation  method  to  pre- 
dict detachment,  even  though  the  prescribed  pressure  gradient  was  ob- 
tained experimentally  from  a separated  flow  in  which  pressure  gradient 
relaxation  has  occurred.  This,  too,  may  be  observed  in  several  of  the 
predictions  in  [3],  and  is  again  the  result  of  not  including  the  free- 
stream  interaction  explicitly  into  the  calculation. 
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Fig.  D.  Behavior  of  6 in  cyclic  iteration. 

Because  of  the  unavoidable  approximations  involved  in  modeling  the 
turbulent  shear  stresses  and  small  errors  in  measurement,  the  phase  re- 
lationship between  the  pressure  gradient  and  the  dependent  b.l.  variables 
in  the  actual  flow  can  never  be  exactly  duplicated  in  a calculation  using 

this  very  same  pressure  gradient  as  a boundary  condition.  Therefore  the 

* 

experimental  growth  of  6 does  not  exactly  match  the  calculated  value. 

If  the  calculated  value  is  slightly  ahead  of  the  measured  one,  runaway 
* * 
growth  of  6 will  occur.  Conversely,  if  the  calculated  6 lags, 

the  freestream  pressure  gradient  will  relax  prematurely  and  the  b.l.  will 

not  detach;  near  detachment,  the  classical  b.l.  procedure  tends  to  become 

unstable . 

Several  methods  are  used  in  practice  to  avoid  this  problem  of  non- 
detachment of  a calculated  b.l.  from  experimental  data.  A very  popular 
scheme  is  the  "frozen  dp/dx"  method,  wherein  the  pressure  gradient  is 
maintained  at  its  maximum  value  and  prescribed  on  the  b.l.  until  it  de- 
taches. Cebeci  et  al . |6|  present  several  comparisons  of  this  method 
against  the  separation  criteria  of  Head,  Goldschmied  and  Stratford. 
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There  are  several  objections  to  the  use  of  methods  such  as  the  fro- 
zen dp/dx  method.  A primary  one  is  that  it  will  predict  detachment  in 
cases  where  none  should  occur,  such  as  in  the  case  of  a flow  which  is 
decelerated  and  then  allowed  to  relax.  In  addition,  there  is  no  physical 
basis  for  the  method,  even  though  it  gives  good  answers  for  the  location 
of  zero  wall  shear  for  rapidly  detaching  flows. 


C.  Simultaneous  Iteration 

The  heuristic  explanation  given  above  suggests  that  the  "separation 
singularity"  and  the  inability  to  predict  detachment  is  nothing  more  than 
prescription  of  the  wrong  boundary  conditions  on  the  b.l.  equations. 

If  a method  could  be  devised  wherein  the  pressure  gradient  at  any 
given  point  is  the  result  of  jiutual  interaction  between  the  b.l.  and  the 
inviscid  core,  then  no  such  singular  behavior  should  occur.  In  this  type 
of  calculation,  the  pressure  gradient  (or  equivalently  the  core  velocity 
at  the  edge  of  the  b.l.  , u ) is  assumed  unknown,  and  an  additional  equa- 
tion, commonly  a 1-D  continuity  equation  in  the  core,  is  added.  This  set 
of  equations  is  solved  simultaneously  at  each  step  along  the  flow.  We 
shall  call  this  scheme  "SIMULTANEOUS  ITERATION";  its  main  features  are 
outlined  in  Fig.  E.  In  Chapter  Five  the  method  will  be  extended  to  the 
case  where  the  edge  velocity  is  obtained  from  a solution  of  the  2-D  La- 
Place's  equation  in  the  inviscid  core. 

A mathematical  description  of  cyclic  and  simultaneous  iteration 
follows.  For  steady,  two-dimensional,  incompressible  flow,  the  b.l. 
equations  are 


x momentum: 


3u  jhj  _ _1  3j)  _3 t_ 

3x  v 3y  P 3x  3y 


(1-1) 


continuity : 


y momentum: 


l£  * 
3y 


o 


3u  Sv_ 
3x  3y 


0 , 


so  that 


i i_E. 
p 3x 


u 

oo 


du 


oo 


dx 


(1-2) 

(1-3) 


Consider  cyclic  iteration  number  n.  The  dependent  variables  are 
calculated  using  the  pressure  gradient  obtained  in  iteration  (n-1). 
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(1-4) 
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du 

oc 

dx 


(n-l ) 


(n) 


If  the  pressure  gradient  is  approximately  the  same  for  both  itera- 
tions, such  as  for  attached  flows,  then  this  set  of  equations  gives  a 
good  solution.  If,  however,  dp/dx  varies  greatly  between  iterations, 
i.e. , 


(n) 


¥ [pj 


(n-l) 


then  either  one  obtains  the  solution  to  the  wrong  problem,  or  the  equa- 
tions diverge,  giving  no  solution  at  all. 

In  simultaneous  iteration,  the  pressure  gradient  is  replaced  by  that 
at  the  current  iteration,  so  that  all  quantities  are  now  at  step  n. 

That  is,  Eqn.  (1-4)  becomes 


[■ 


9u 

3x 


+ v 


3ul  (n) 

3yJ 


du 


oo 


dx 


+ 


3t 

3y 


(n) 


(1-5) 


[u  du^/dx] ^ is  now  unknown  and  must  be  supplied  from  an  addi- 
tional relationship  involving  the  potential  core.  In  Chapter  Two  a one- 
dimensional core  equation  is  used,  while  in  Chapter  Five  the  method  will 
be  extended  to  include  a two-dimensional  core.  In  either  case,  all  quan- 
tities are  expressed  in  terms  of  values  at  step  n.  In  effect,  we  have 
converted  from  an  explicit  to  an  implicit  set  of  equations,  with  a cor- 
responding gain  in  numerical  stability. 

A clear  explanation  of  simultaneous  iteration  as  applied  to  a dissi- 
pation integral  type  b.l.  calculation  can  be  found  in  Gerhart  [50]. 


D.  Previous  Work 

Moses  [7]  used  the  simultaneous  iteration  concept  to  calculate  the 
flow  in  incompressible  diffusers  with  detached  b.l.'s.  Even  though  he 
used  a power-law  velocity  profile  and  rather  simple  prediction  schemes, 
he  was  able  to  obtain  fair  agreement  with  data  for  diffusers  operating 
in  the  early  portions  of  the  transitory  stall  regime.  The  most  signifi- 
cant contribution  of  his  work  was  to  recognize  the  need  for  including  a 
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1-D  core  equation,  thereby  enabling  him  to  avoid  the  singular  behavior 
of  the  equations  near  detachment. 

Moses  was  much  criticized  for  having  the  audacity  to  attempt  the 
calculation  of  flows  at  and  beyond  detachment  [8].  Unfortunately,  some 
of  the  correlations  used  by  him  were  questionable  in  the  light  of  the 
then  available  data.  The  discussion  of  his  work  in  the  literature  quickly 
became  bogged  down  in  arguments  over  the  validity  of  details  of  these 
correlations,  and  the  central  idea,  that  of  simultaneous  iteration,  was 
largely  forgotten. 

The  next  few  years  saw  tremendous  activity  in  the  development  of 
prediction  methods  for  turbulent  b.l.'s,  as  exemplified  by  the  1968  Stan- 
ford Conference  [3]  on  Computation  of  Turbulent  Boundary  Layers.  In  all, 
28  methods  for  the  calculation  of  TBL's  with  prescribed  pressure  gradient 
were  presented.  These  varied  from  simple  correlative  integral  methods  to 
rather  complicated  differential  methods  using  sophisticated  turbulence 
closure  schemes.  None  of  the  methods  presented  was  able  to  calculate 
separating  flows  near  detachment  adequately.  This  is  not  surprising, 
since  they  were  all  calculated  with  prescribed  pressure  gradient  without 
taking  the  freestream  interaction  into  account. 

It  is  apparent  from  the  discussions  that  some  predictors  were 
acutely  aware  of  the  need  for  including  this  interaction  for  separating 
flows.  Nevertheless,  the  majority  of  attendees  bypassed  this  in  favor 
of  discussions  involving  the  validity  of  the  b.l.  equations,  the  contri- 
butions from  normal  stresses,  curvature  effects,  etc. 

The  controversy  regarding  the  inability  of  the  b.l.  equations  with 
prescribed  pressure  gradient  to  predict  detaching  flows  is  still  raging. 

As  late  as  1975,  attendees  at  the  AGARD  Separating  Flow  Conference  [9] 
were  still  debating  the  same  questions  as  at  the  '68  Stanford  Conference. 
The  idea  of  simultaneous  iteration  being  the  key  to  removing  the  singular 
behavior  near  detachment  is  still  far  from  being  universally  accepted. 

In  1972,  Bower  [10]  extended  Moses'  calculation  to  include  compres- 
sible flow  in  axisymmetric  diffusers.  He  retained  the  1-D  core  assump- 
tion and  used  a dissipation  integral  b.l.  method.  The  dissipation  inte- 
gral, 0^  was  related  to  the  shape  factor  H through  an  empirical 
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correlation  due  to  Alber  [11].  The  energy  shape  factor,  H,  w.  elated 
to  H through  the  Escudier-Nicoll  correlation  [3], 


H = 1.431  - .0971/H  + .775/H 


(1-6) 


A limiting  form  of  Coles'  velocity  profile  was  used,  with  Re.  -*■  <» 

o 

giving  a one-parameter  family. 

£ = + *n  (f)An  (-565  + 1 (nrH1  - cos  11  si 

d-7) 

Skin  friction,  C^/2,  was  obtained  through  the  Ludweig-Tillmann  corre- 
lation. 


= 0.123  Re”’268  10-678  H 


(1-8) 


Bower's  predictions  for  the  diffusers  operating  in  the  early  por- 
tions of  the  transitory  stall  regime  are  quite  good.  Nevertheless,  his 
calculation  method  can  be  criticized  on  several  grounds. 

The  one-parameter  velocity  profile,  Eqn.  (1-7),  is  a poor  represen- 
tation of  actual  b.l.'s  in  adverse  pressure  gradients,  even  though  it 
does  permit  backflow.  The  empirical  H vs.  H relationship,  Eqn.  (1-6), 
is  valid  only  for  1.25  £ H <_  2.8  (Ref.  [3],  pp.  136-138),  but  is  used 
in  this  method  for  H up  to  12.0.  The  Ludweig-Tillman  correlation, 

Eqn.  (1-8),  is  always  positive,  so  that  zero  or  negative  wall  shear  val- 
ues cannot  be  represented,  no  matter  how  large  the  values  of  Re.  and 
H.  As  a result,  the  location  of  zero  wall  shear  cannot  be  determined, 
and  the  rather  arbitrary  value  of  H = 1.8  was  used  as  an  indication  of 
detachment . 

However,  it  is  well  known  that  detachment  is  not  a unique  function 

* 

of  H,  being  in  fact  a stronger  function  of  the  blockage,  6 . This  is 
apparent  in  the  work  of  Sandborn  and  Kline  [51],  who  postulate  the  begin- 
ning of  intermittent  detachment  at  a point  where  H = H , where 


1 -6*/6 


(1-9) 


11 


I 


.iZTif  3m. 


Also,  Senoo  and  Nishl  [13]  obtained  an  empirical  stall  limit  rela- 
tion for  diffusers. 


I 


Hg  = 1.8  + 3.75  B , (1-10) 

* 

where  B is  the  local  value  of  the  blockage  factor,  26  /W. 

The  use  of  any  empirical  correlation  to  determine  detachment  is 
clearly  undesirable,  since  it  limits  the  probable  generality  of  the  pro- 
cedure and  is  hence  to  be  avoided  if  possible. 

In  view  of  these  residual  difficulties  in  the  work  of  Bower,  the 
relatively  good  results  achieved  strongly  support  (but  do  not  definitely 
prove)  the  idea  that  the  central  difficulties  in  predicting  detachment 
and  separated  flows  can  be  cured  or  strongly  alleviated  by  simultaneous 
iteration.  To  put  this  differently,  as  already  found  by  Woolley  [4]  and 
White  [5]  for  fully  stalled  flows,  the  crucial  matter  is  to  get  the  in- 
teraction between  the  blockage  effects  of  the  separated  zone  and  the 
outer  flow  modeled  adequately;  all  other  effects  are  less  important  to 
adequate  predictions.  What  is  suggested  here,  then,  is  that  the  same  is 
true  of  detachment  and  detaching  flows,  and  that  for  such  cases  simulta- 
neous rather  than  cyclic  iteration  is  necessary.  It  is  this  idea  plus 
the  specific  details  needed  to  alleviate  the  problems  relative  to  Bower's 
work  that  are  central  to  the  work  that  follows. 
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CHAPTER  TWO 

UNIFIED  INTEGRAL  METHOD 
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A.  Requirements  for  Boundary  Layer  Prediction  Method 

The  survey  of  currently  available  prediction  methods  for  turbulent 
b.l.'s  presented  in  the  last  chapter  showed  that  calculation  methods  for 
attached  b.l.'s  are  highly  developed.  The  converse  is  true  for  detached 
flows  and  b.l.'s  that  are  in  the  process  of  detaching,  both  of  which  must 
be  calculated  in  simultaneous  fashion  with  the  bounding  freestream.  It 
was  felt  that  a new  calculation  method  was  needed  to  be  able  to  extend 
the  diffuser  calculations  deeper  into  the  transitory  stall  regime.  The 
requirements  for  such  a method  are  that: 

(a)  The  equations  simultaneously  solve  for  the  boundary  layer  and 
the  freestream. 

(b)  The  velocity  profile  family  be  capable  of  representing  both 
attached  and  detached  flows. 

(c)  The  auxiliary  equation,  turbulent  shear  stress  model  and  its 
associated  correlations  be  valid  for  attached  and  detached  flows. 

(d)  The  set  of  equations  should  not  introduce  any  singularities  at 
the  detachment  point. 

(e)  Detachment  should  occur  "naturally"  and  be  perceptible  as  hav- 
ing occurred  without  recourse  to  any  empiricism  such  as  the  frozen 
dp/dx  method  or  a detachment  criterion.  That  is,  the  desirable  detach- 
ment criterion  is  = 0 (on  the  average). 

(f)  The  method  should  be  fast,  since  we  expect  to  use  it  in  an 
iterative  fashion. 

(g)  The  core  velocity  should  be  obtained  from  a solution  of  the 
elliptic  irrotational  core,  so  as  to  include  downstream  effects  on  the 
upstream  flow. 

The  rest  of  this  chapter  develops  such  a method,  the  "Unified  Inte- 
gral Method  " (UIM) , with  a 1-D  core.  The  extension  to  the  2-D  case  is 
deferred  to  Chapter  Five. 
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B . Integral  Methods  in  General 

A brief  summary  of  integral  b.l.  methods  will  be  presented  before 
proceeding  with  the  development  of  the  UIM  equations. 

All  integral  methods  use  the  von  Kdrman  momentum  integral  equation, 

d0  4.  ,9, m 6 dU°°  _ Cf  . 1 f6  3 (o  n 

d x+(2+H)^17  “ T + ~2  J 9^(u  "v  )dy  > (2_1) 

00  u * o ' ' 

co 

The  normal  stress  term  is  usually  neglected,  although  there  is  some 
evidence  that  its  value  may  be  large  near  detachment.  This  equation  may 
be  parametrically  expressed  as 


= f ]_(H»  e,  u.,  Cf/2) 


(2-2) 


Consider  the  case  of  a prescribed  pressure  gradient  calculation 
where  u is  a known  function  of  the  streamwise  coordinate  x.  Two 

OO 

more  equations  are  needed  to  solve  for  the  three  unknowns  0,  H,  and 
C^/2.  The  differences  in  the  various  methods  arise  in  the  procedure 
used  to  close  the  set  of  equations. 

Many  methods  use  an  empirical  equation  relating  the  skin  friction 
C^/2  to  the  calculation  variables.  The  most  commonly  used  is  the 
Ludweig-Tillmann  correlation. 


= 0.123  Re'*128  10-678  H 

L o 


(2-3) 


The  last  equation  remaining  is  called  the  "auxiliary"  equation  and 
relates  the  growth  of  the  shape  factor  H to  the  other  b.l.  parameters. 
One  method  of  obtaining  this  equation  is  by  taking  moments  in  u or  y 
of  the  two-dimensional  b.l.  equations  before  integrating  across  the  layer. 
The  first  moment  in  u gives  the  "mechanical  energy"  equation. 


6 4^  = (H-l)  ~~  — H ~7t  + C , 

dx  u dx  2 D 


( 2-4) 


where 


2 f du  . 

~ij  1 ^ dy 
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(2-5) 
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The  first  moment  in  y gives  the  "moment  of  momentum"  equation, 


r6  f 3u2  3 ( fy  5u  , \]  , fi2 

J0  uJo  ^dyl  dy  = t u“ 


(2-6) 


where 


(2-7) 


Additional  unknowns  H,  C^,  CT  have  appeared  in  both  auxiliary  mo- 
ment equations  (2-4)  and  (2-6),  and  these  must  be  related  back  to  the 
primary  variables  H,  9,  and  C^/2.  At  this  stage  a model  equation  for 
the  turbulent  shear  stresses  and  a velocity  profile  family  must  be  intro- 
duced. The  turbulence  model  relates  CL  or  C to  the  mean  flow  param- 

D T _ 

eters,  while  the  velocity  profile  family  allows  H to  be  expressed  in 
terms  of  H for  Eqn.  (2-4)  and  permits  Eqn.  (2-6)  to  be  integrated.  For 
details  of  this  process,  see  the  review  papers  by  Reynolds  [3],  Rotta 
[14],  and  the  introductory  sections  of  Hirst  and  Reynolds  [15]. 

Head  [16]  used  the  growth  rate  of  the  turbulent-nonturbulent  front 
to  derive  an  auxiliary  equation.  The  rate  at  which  the  b.l.  spreads  into 
the  irrotational  fluid  is  the  entrainment  rate  dQ/dx  and  may  be  ex- 
pressed as  a function  of  a new  shape  factor  Hg_^*  , 


-/• 


dy  = u (6-6  ) 


(2-8) 


W<5*’  u»’  6-6  } 


(2-9) 


where 


H6-6* 


(2-10) 


H.  r*  is  in  turn  related  back  to  H through  another  correlation, 

o-o 

closing  this  set  of  equations. 

A survey  of  the  literature  will  show  the  large  variety  of  auxiliary 
equations  that  have  been  used.  This  is  a consequence  of  the  fact  that 
no  "exact"  independent  equation  is  available.  It  is  therefore  important 


to  understand  exactly  what  the  auxiliary  equation  provides  in  the  way  of 
new  information. 

We  note  that  there  is  no  term  involving  the  turbulent  shear  stresses 
in  the  momentum  integral  equation  (2-1).  Therefore,  the  most  important 
function  of  the  auxiliary  equation  is  to  supply  information  regarding  the 
shear  stresses  in  the  b.l.  The  second  requirement  is  that  it  truly  con- 
tain independent  information.  For  instance.  Hirst  et  al.  [15]  found 
that  the  mechanical  energy  equation  may  not  be  completely  independent  of 
the  momentum  integral  equation.  This  is  so  since  u is  fairly  constant 
across  the  layer,  and  the  resulting  set  of  equations  is  almost  redundant. 

Studies  by  Hirst  et  al.  [15]  and  Thompson  [35]  showed  that  the  en- 
trainment method  of  Head  appeared  to  work  better  than  other  available 
methods  for  a large  variety  of  b.l.'s.  They  hypothesized  that  perhaps 
this  technique  contained  "more"  independent  information  regarding  the 
turbulence.  We  shall  therefore  use  the  entrainment  concept,  but  extend 
its  applicability  to  enable  calculation  of  detached  flows. 

To  summarize,  the  auxiliary  equation  is  of  the  form 

^ = f 2 (H , 0,  Uoo,  Cf / 2 , <]>(T))  . (2-11) 

The  turbulence  model  equation  is  of  the  form 

<Kt)  = f3(H,  0,  uOT,  Cf/2)  . (2-12) 

<{>(t)  is  a functional  representation  of  shear  stress  integrals  such 
as  0^  or  C^.  in  Eqns . (2-4)  and  (2-6).  The  closure  model  for  cf) ( T ) 
relates  it  back  to  known  quantities  such  as  H,  0,  u^,  and  C^/2,  as 
shown  through  Eqn.  (2-12). 

The  skin  friction  equation  is  obtained  from  a correlation  of  the 

form 

= f4(0,  H,  uj  . (2-13) 

through  (2-13)  permit  the  b.l.  parameters 


to  be  calculated  in  a stepwise  marching  fashion  along  the  flow. 


Cf/2 


Equations  (2-2)  and  (2-11) 


As  mentioned  before,  all  of  the  above  methods  work  quite  well  for 
accelerating  and  flat-plate  flows,  and  reasonably  well  for  decelerating 
flows  which  are  far  from  detachment.  Neither  the  velocity  profile  fam- 
ily nor  the  auxiliary  equations  are  valid  at  or  beyond  detachment.  We 
proceed  therefore  to  tailor  the  UIM  to  be  able  to  do  this  by  first  exam- 
ining a velocity  profile  family  that  is  capable  of  representing  both 
attached  and  detached  flows,  and  then  developing  an  auxiliary  equation 
that  works  over  this  entire  range. 


C.  Velocity  Profile  Family 

It  is  generally  accepted  that  typical  TBL  velocity  profiles  can  be 
represented  by  the  combination  of  an  inner-wall-dominated  layer  plus  an 
outer  "wake-like"  structure.  One  such  velocity  profile  family  is  the  log 
law  of  the  wall  matched  to  Coles'  [17]  "law  of  the  wake"  outer  profile. 


(2-14) 


where 

k is  the  von  Karman  constant  = 0.41, 
c is  the  wall  constant  = 5.0, 
u is  the  shear  velocity  = /t  /p , 

T J CO 

n is  the  wake  amplitude. 


This  profile  gives  excellent  results  for  attached  flows,  but  cannot 
be  used  without  modification  for  either  detached  or  detaching  flows.  The 
difficulty  in  representing  detaching  flows  may  be  seen  by  taking  the  limit 
as  u^  -*■  0 after  setting  u = u^  at  y = 6. 


i.e. , 


u 


00 


— lim  (nu  ) 
k _ T 
u -*-0 
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So  II  ->■  “ as  u ->  0 in  order  to  keep  the  limit  finite. 

T 

Beyond  detachment,  the  wall  shear  x^  is  negative  and  is  not 

even  defined. 

However,  a simple  modification  of  (2-14)  together  with  redefinition 

of  u for  reversed  flows  permit  the  desired  representation, 
x 

Let 


and 


ut  £ (sgn 


+ A y'V 


The  modified  velocity  profile  family  is 


y |u 


+ c 


+ -1  (l  - cos 
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(2-15) 


where  uQ  is  the  redefined  wake  amplitude,  and  c is  a new  constant, 

p 


c = ck  = 2.05. 
Define 


A 
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A 
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(2-16) 


n = 


£ X 
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At  the  edge  of  the  b.l.,  n = 1,  and  u = u . 
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(2-17) 


Subtracting  Eqn.  (2-17)  from  (2-15)  and  substituting  from  Eqn. 
(2-16),  we  get  the  desired  profile, 
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— = 1 + V In  n - V cos2 
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(2-18) 


Fig.  F shows  the  contribution  of  each  of  the  terms  in  this  equation 
to  the  velocity  profile. 

We  note  from  Eqn.  (2-18)  and  the  above  sketch  that  uQ  is  normally 

p 

positive  and  that  u„  < u for  attached  flows  while  u„  > u for  de- 

v 8 » 3 oo 

tached  flows.  This  form  of  the  velocity  profile  has  been  used  by  McDon- 
ald and  Stoddard  [18]  and  Nash  and  Hicks  [3]  for  attached  flows.  Kuhn 
and  Nielsen  [19]  attempted  to  calculate  detached  flows  using  this  pro- 
file. The  ability  of  Eqn.  (2-15)  to  represent  attached  and  detached 
flows  is  shown  in  Figures  2 and  3. 

Alber  et  al.  [20]  extended  the  applicability  of  this  profile  to 
represent  compressible  flows,  and  concluded  that  the  Coles  type  formula- 
tion is  an  adequate  representation  of  the  velocity  field  both  upstream 
and  downstream  of  detachment.  The  measured  detachment  profile  does  not, 
however,  correspond  to  Coles'  zero  wall-friction  profile.  Over  most  of 
the  flow,  however,  a fair  to  good  fit  with  experimental  data  was  obtained. 
It  should  therefore  be  adequate  for  use  in  an  integral  prediction  method. 

Using  the  velocity  profile,  Eqn.  (2-18),  it  is  possible  to  develop 

* 

relationships  among  the  b.l.  integral  parameters,  6 , 0,  and  V , V . 

* IB 

Substituting  Eqn.  (2-18)  into  the  definitions  of  6 and  0 and  on 
integrating  across  the  b.l.,  we  get 


* V 

i_  « v 

6 T 2 


= 2V2  + | V2  + 1.58949  VTVB 


(2-19) 


(2-20) 


These  two  equations  give  an  unambiguous  definition  of  6. 


Momentum  Integral  Equation 

The  momentum  integral  equation  for  steady,  2-D,  incompressible  flow 

f + <2+“>  f TT  ' T * -r/  T*  (u'2-'',2>  ■ <2-21) 

on  n ^ n 
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' * 


This  may  be  rewritten  in  terms  of  the  proposed  dependent  variables 


6,  u0,  u , and  U by  defining 

p T 00 


/ £ 

N u2J 

U * o 


3 < 1 2 .2%  . 

— (u1  -v'  ) dy 
dx 


(2-22) 


and  noting  that 


* - G)  • 


(2-23) 


Differentiating  (2-19)  with  respect  to  x, 

* * . du  . / u uQ  \ dU  k du 

d5  6 d6  _6 t _6_  /_L  + _JL  I ” + A §.  rj.ou\ 

dx  ' 6 dx  + <Vx  dx  " Um  \kUco  2uJ  dx  2U<jo  dx  * (2_  > 


Similarly,  differentiating  (2-20)  and  rearranging  gives 


de  = 1 1 

dx  \ 


2 3 2 

/ _ — v 

T 8 B 


2 

. /4Y  . 3 „2  6 . 3.17898 

+ 1— + 4 VB  r + —3-  u 
\ 00  00  <u 

00 

/ 6 3 6u8  1.58949  r\ 


V \ 

/46 

' B/  dx 

f.  _ 

6vt 

eV 

u 

00 

dUg 

dx 

„ 1.58949V_ 6 , 

/2  + L - -1 

T kU  kL 


J dX 


(2-25) 


Substituting  (2-22)  through  (2-25)  into  (2-21),  rearranging,  and 
using  (2-19)  and  (2-20)  gives 


(!)idi|  + (f)(MvB-  i-589i9vT)}^f| 

♦ (£)  (*  - «t  - l-58949  '»){£}♦(£){£] 


(2-2(>) 


2 2 

‘ VT  + % 


Equation  (2-26)  is  the  form  of  the  momentum  integral  equation  used 
in  the  computations.  The  normal  stress  term,  V^,  has  been  carried 
along  for  completeness  and  is  not  used  further  in  this  investigation. 
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1.  Outer-Edge  Matching  Equation 


Differentiating  Eqn.  (2-17)  in  the  streamwise  direction  and  rearrang- 


ing gives 


■H 


' 6 lUr I \ ( dur)  (SSn  UT)  / d I UT  I 

*n-7~  + c ) Ur  + “ ^ 6^r+ 


•g  £)♦(£)  • 


(2-27) 


(sgn  u ) I u^ | 


d u I du 

T T 

(sgn  u ) — z = -3 — . 

0 x dx  dx 


(2-28) 


(2-29) 


Equation  (2-29)  is  valid  everywhere  except  in  the  case  where  u^ 
changes  sign.  There  is  then  an  ambiguity  in  the  sign  of  the  resulting 
value  of  u^,  which  may  be  resolved  by  using  physical  insight  from  the 
velocity  profile.  When  u^  > u^,  then  u^  < 0 and  vice-versa. 

Therefore, 


I uT  | sgn(ug-Uoo) 


(2-30) 


Rearranging  Eqn.  (2-17),  we  get 


1 / *l\l  .\  U«=  - 

- I Hn  — + c I = 

K \ 1 


Substituting  (2-28)  through  (2-31)  in  (2-27)  gives 


id«i  . , . idu6i  . r 

!dx  j + (ux}  1171  + 


(v  + u»  - ub)  ["57} + (_ut)  (“srj 


(2-31) 


(2-32) 


F.  Entrainment  Equation 

The  concept  of  entrainment  will  be  used  to  derive  the  auxiliary 
equation.  The  calculation  method  of  Bradshaw  et  al.  [21]  uses  a corre- 
lation between  the  nond  imens Iona  1 entrainment  rate 

1 d * 

TT  f 0 (A-A  ) 

U d x * 


22 


and  the  maximum  shear  stress  in  the  b.l.. 


t /pi!2  . 
max  <» 


This  correlation 


has  been  revised  in  Fig.  A to  include  data  from  recent  experiments,  and 
shows  that  the  entrainment  rate  is  about  ten  times  the  maximum  shear 


stress.  That  is, 


-L  Jl|  „ 

U dx  I “ 

oo  1 — 


(6-6*) 


10  t /pU 
max  “ 


(2-33) 


The  remarkable  feature  of  this  correlation  is  that  it  seems  to  apply 
to  both  attached  and  detached  flows.  It  works  equally  well  for  b.l.'s  in 
favorable  or  adverse  pressure  gradients,  provided  that  for  accelerating 
flows  the  maximum  shear  stress  is  evaluated  at  n = y/6  = 1/A,  despite 
the  fact  that  the  maximum  shear  stress  for  an  accelerating  b.l.  occurs 
at  the  wall. 

Differentiating  Eqn.  (2-33)  in  the  x direction  and  substituting 
* 

for  d<5  /dx  from  (2-2A)  gives 


= 10  T /pU 

max  00 

(2-3A) 


We  require  x / p to  be  able  to  use  (2-3A).  The  distance  from  the 

max 

wall  y/6  at  which  the  shear  stress  is  maximum  will  be  obtained  from 
another  correlation.  The  velocity  profile  can  be  differentiated  and 
evaluated  at  this  y/6  location  to  give  the  value  of  3u/3y  correspond- 
ing to  maximum  x.  It  is  then  possible  to  compute  x /p  by  using  an 

ms  x 

eddy-viscosity  model. 

A plot  of  (u/U  ) at  which  the  maximum  shear  stress  occurs 

1 OO  x 


(Fig.  G)  as  a function  of  the  ratio  of  the  wall  to  wake  velocities, 

2u  /ku  , is  shown  in  Fig.  5.  There  is  a fair  amount  of  scatter  and  a 
x p 

clearer  picture  emerges  when  only  equilibrium  b.l.'s  and  detached  flows 
are  plotted  (Fig.  6).  It  is  apparent  that  the  velocity  ratio  at  which 
x occurs,  denoted  by  (tt*)  , may  be  quite  well  represented  by 


equilibrium  b.l.'s: 


(r)  • 

\ oo  f J 

max 

(f), 

V oo  / m 


= 0.76 


detached  flows: 


(il 


00 ' max 
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f(nT  ) - VT  in  n, 

' max 


,,  2 u , , 

- VB  COS  2 \ + 1 

max  max 


" (Ot. 


(2-35) 


This  gives  n • Differentiating  Eqn.  (2-18),  we  get 
Tmax 


= ^ /— 
\3y/T  6 \ T) 

' * ' lmax  v x 


U“/  ^ + V f 

Tin + — Sln  *n 


t ) 


(2-36) 


On  substituting  n_  we  obtain  (3u/3y)T  . This  may  now  be 

max  max 

used  through  an  eddy  viscosity  formulation  to  obtain  the  maximum  shear 

stress  t ; i.e.,  we  assume 

max 


3y  ’ 


(2-37) 


where  e is  an  eddy  viscosity. 

For  the  outer  portion  of  equilibrium  b.l.'s  (y/6  > 0.2),  Clauser 
[22]  showed  that  e may  be  approximately  represented  by 


K U 6 , 


(2-38) 


where  K = .0168. 
e 

For  nonequilibrium  b.l.'s,  the  relationship  is  no  longer  this  simple, 
and  many  efforts  have  been  made  to  obtain  a universal  formulation  for  e . 
McD.  Galbraith  and  Head  [23]  present  an  extensive  summary  of  many  of 
these  attempts  and  compare  the  results  with  experiments. 

Kuhn  and  Nielsen  [19]  included  the  effect  of  pressure  gradient  and 
intermittency  and  obtained 


- - K yU  6*  (f51)  , 

o e o°  \3y/ 


(2-39) 


where 


.013  + .0038  e 


-8/15 


6 - 

t dx 
0) 


and  intermittency  y - (l+9n  ) 
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The  pressure  gradient  parameter  6 is  small  for  mild  pressure  gra- 
dients and  flows  far  from  detachment.  This  allows  to  approach  Clau- 

ser's  value  of  .0168  for  attached  equilibrium  b.l.'s  and  the  limiting 
value  of  .013  for  free  shear  flows  such  as  the  flow  at  a free  jet  bound- 
ary (Schlichting  [2],  pp.  681-707)  for  which  8 -*  00 . For  accelerating 
flows,  8 is  set  equal  to  zero. 

The  maximum  shear  stress  in  an  equilibrium  b.l.  can  thus  be  obtained 


from 


(f). 


= (.013  + .0038 


e-6/15)(l 


SnV1 


U 6 


max , eq 


(*?), 

' ' max 


(2-40) 


For  a nonequilibrium  b.l.,  this  expression  has  to  be  modified  to 
take  into  account  upstream  history  effects.  The  fluid  near  the  wall  in 
a TBL  is  in  local  equilibrium  in  the  sense  that  it  adjusts  very  rapidly 
to  external  changes,  such  as  the  pressure  gradient.  The  outer  layers, 
however,  are  dominated  by  large  eddies  that  have  considerable  inertia, 
so  that  it  has  finite  adjustment  time.  The  outer  layer  therefore  "lags" 
behind  the  local  pressure  gradient.  Typical  behavior  of  the  velocity 
profile  in  response  to  sudden  removal  of  the  pressure  gradient  is  shown 
schematically  in  Fig.  H,  adapted  from  White  [25]. 


PRESSURE 
GRADIENT 
REMOVED 
AT  x = X0 


X0+8 


EQUILIBRIUM 


X0+28 

^ X0+58 

LOGARITHMIC 
OVERLAP  LAYER 


LOG  (y+) 


Fig.  H.  Relaxation  following  removal  of  pressure  gradient. 
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For  the  flow  shown  in  Fig.  H,  the  velocity  profile  takes  56  after 
removal  of  the  pressure  gradient  before  it  reaches  equilibrium.  Another 
example  of  shear  stress  lag  may  be  seen  in  the  relaxing  flow  of  Goldberg 
[26],  as  presented  by  McDonald  et  al.  [18].  Fig.  I shows  the  lag  between 
the  measured  shear  stress  integral  and  its  equilibrium  value  . 

The  calculation  of  shear  stresses  from  an  equilibrium  condition  will 
therefore  give  erroneous  results. 

One  method  of  accounting  for  the  departure  from  equilibrium  is  to 
relate  the  equilibrium  and  nonequilibrium  values  through  a first-order 
differential  equation,  commonly  called  a lag  equation. 


d / Tm ax  \ _ A_  /J max,  eg  _ Tmax  \ 
dx  \ p / 6\p  P / ‘ 


(2-41) 


The  lag  parameter  A is  obtained  from  numerical  experiments. 

It  should  be  noted  that  the  lag  equation  does  not  model  a primary 
term,  but  only  corrects  a deviation  of  what  would  otherwise  be  an  error 
in  a primary  term.  Since  this  deviation  is  usually  small  and  only  sig- 
nificant for  rapid  changes  in  the  "environment"  of  the  shear  layer  in 
the  streamwise  direction,  the  form  of  the  lag  equation  used  is  not  criti- 
cal. Hence  a simple  first-order  diffusion  equation  should  be  sufficient. 
That  this  is  so  is  demonstrated  by  the  results  in  the  1968  Conference 
[3]. 

A summary  of  the  process  used  to  obtain  the  right-hand  side  of  the 
entrainment  equation,  which  is  now  expressed  entirely  in  terms  of  known 
quantities,  is  shown  in  Fig.  7.  Fig.  8 compares  Tmax^P  measured  by 
Strickland  and  Simpson  [32]  with  that  from  Eqn.  (2-40).  The  measured 
entrainment  rate  is  also  shown.  The  agreement  is  quite  good  for  the 
attached  flow,  and  the  last  few  detached  flow  points,  and  fair  to  poor 
for  the  rest.  Except  for  the  last  station,  Tmax/p  *s  overpredicted . 
This  is  consistent  with  our  expectations,  since  these  values  were  calcu- 
lated from  the  measured  mean  velocity  profile  and  correspond  to  the 
equilibrium  case.  Shear  lag  will  decrease  these  computed  values.  The 
worst  match  is  at  station  124.3,  the  location  of  intermittent  detach- 
ment, where  the  uncertainty  in  both  the  measured  and  calculated  values 
is  the  greatest. 
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Fig.  I.  Behavior  of  equilibrium  shear  stress  integral  C and 

measured  value  CT  for  the  relaxing  flow  of  Goldberg  [26]. 
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To  calculate  a b.l.  with  prescribed  pressure  gradient,  terms  involv- 
ing dU^/dx  are  moved  to  the  right-hand  side  of  Eqns.  (2-26),  (2-32), 
and  (2-34),  which  may  then  be  integrated  in  a stepwise  fashion  along  the 
flow. 

For  detaching  flows  that  must  be  calculated  by  simultaneous  itera- 
tion, an  additional  equation  is  needed  since  is  now  an  unknown.  In 

this  chapter  we  shall  use  a 1-D  continuity  equation  across  the  diffuser 
width. 

G.  One-Dimensional  Core  Equation 

Consider  flow  in  a diffuser  of  width  W(x)  with  a uniform  1-D 
velocity  distribution  in  the  core  (Fig.  J) . 


8* 


8< 


Fig.  J.  1-D  core  diffuser  nomenclature. 

Assuming  symmetrical  b.l.'s  and  ignoring  end-wall  effects,  the  con- 
tinuity equation  at  any  section  x is 


U (W  - 26  /cos  0)  . 


(2-42) 
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On  differentiating  (2-42)  in  the  x direction,  substituting  for 


d6  /dx  from  (2-24),  and  manipulating,  we  get 


/-6*\  I d6  ) / -6  \ (du6  ) / -i 5 \ (dUT  ) [“cos  0 26*  \ 6* 

(— ) tel + (ir)  + (^r)  j^rj + ~ ^tt)  + r 


dU 

oo 

dx 


cos  6 dW 
2 dx 


(2-43) 


H . Solution  of  the  UIM  Equations 

The  addition  of  Eqn.  (2-43)  to  (2-26),  (2-32),  and  (2-34)  closes  the 
set.  These  may  now  be  written  as  a 4 x 4 matrix  equation  at  each  step 
along  the  flow. 


— — 

1 d6  1 

all  a12  al3  a14 

dx 

dun 

£3  , • • • 

B 

b„ 

21 

dx 

2 

du 

— 

* 

T 

a31 

dx 

dU 

* 

j*41  ’ * a44^ 

dx 

b4 

(2-44) 


The  unknown  x derivatives  are  first  uncoupled  using  Gauss  elimi- 
nation and  the  resulting  set  of  ODE's  solved  with  a fourth-order  Adams- 
Bashforth  predictor-corrector  routine. 


At  detachment,  u^  = 0,  and  both  sides  of  (2-32)  become  identically 


equal  to  zero.  To  prevent  the  coefficient  matrix  from  becoming  singular, 
this  equation  is  removed  for 


u < .025  and  the  reduced  set  solved 
' x ' 


with  the  value  of  du^/dx  frozen  at  the  predetachment  value.  The  re- 


sults are  negligibly  affected  by  varying  this  threshold  value  of 


between  .005  and  .050. 

In  the  next  few  chapters,  the  TBLPM  developed  herein  will  be  used 
to  predict  diffusers  after  being  calibrated  by  comparing  its  performance 
with  that  of  a large  number  of  currently  available  calculation  methods. 
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CHAPTER  THREE 

RESULTS  - BOUNDARY  LAYER  CALCULATIONS  WITH  PRESCRIBED  PRESSURE  GRADIENT 


A.  Determination  of  Lag  Parameter 

The  lag  parameter  X is  still  not  known.  The  relaxation  times  for 
attached  and  detached  flows  are  expected  to  be  quite  different,  since 
the  former  is  wall-dominated  while  the  latter  is  inertia-bound.  Two  lag 
parameters,  X^  for  attached  and  X^  for  detached  flows  were  therefore 
used.  It  is  expected  that  the  "time  constant"  for  detached  flows  will  be 
vefy  short  compared  to  that  for  attached  flows,  since  the  effect  of  the 
wall  is  negligible  and  the  large  fluctuations  present  here  tend  to  rapidly 
destroy  any  upstream  history  effects. 

It  must  be  noted  that  the  existence  of  a lag  between  the  equilibrium 
and  actual  shear  stress  values  for  detached  flows  is  a hypothesis  only, 
and  that  not  enough  data  are  available  to  rule  out  or  confirm  its  exis- 
tence. This  is  in  sharp  contrast  to  the  corresponding  case  for  attached 
flows,  for  which  shear  lag  is  a well-established  phenomenon.  The  results 
of  the  present  investigation  are  also  ambiguous  in  this  regard  and  do 
not  conclusively  support  or  rule  out  the  necessity  for  using  a lag  equa- 
tion for  detached  flows. 

Both  X and  X were  varied  independently  and  the  resulting  pre- 
a d 

dictions  compared  with  data.  The  effect  of  varying  X for  the  relaxing 

SL 

flow  of  Bradshaw  et  al.  [49]  is  shown  in  Figs.  9a  and  9b.  Fig.  20  shows 
the  effect  on  the  unstalled  diffuser  flow  of  Carlson  and  Johnston  [27]. 
Varying  X,  while  keeping  X constant  is  shown  in  Fig.  21  for  a dif- 

Q SL 

fuser  in  transitory  stall,  also  measured  by  the  same  experimenters.  The 
results  of  not  using  any  lag  equation  is  shown  in  all  the  above  cases. 

For  attached  b.l.'s,  varying  X from  .015  to  .035  has  negligible 

d 

effect  on  the  flow.  Not  using  a lag  equation,  however,  does  cause  a 

small  but  discernible  deviation  from  the  data. 

Similar  effects  are  seen  for  detached  flows.  As  X is  varied 

d 

between  0.3  and  0.7,  the  exit  velocity  ratio  u/u_  changes  from  0.68 

Khr 

to  0.64.  This  change  is  of  the  order  of  the  experimental  uncertainty. 

Not  using  any  lag  again  causes  a small  but  detectable  overprediction  of  Cp. 
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We  conclude  that  a lag  equation  is  necessary  for  accurate  predic- 
tions, although  there  is  a large  Latitude  in  the  choice  of  the  numerical 
values  for  \ and  The  numbers  finally  used  were  those  giving  the 

best  match  with  a large  number  of  flows,  and  are 


attached  lag  parameter 


These  results  are  consistent  with  the  remarks  concerning  lag  equa 
tions  following  Eqn.  (2-41). 


Experimental  Data  for  Turbulent  Boundary  Layers 


The  best  collection  of  TBL  data  is  the  extensive  compilation  of 
Coles  and  Hirst  [44] . These  data  have  been  reduced  in  a consistent  man- 
ner; moreover,  the  results  of  a large  number  of  TBLPM's  using  prescribed 
pressure  gradient  to  compute  these  data  are  presented  in  Kline  et  al.  [3 
Any  new  calculation  method  must  be  able  to  predict  satisfactorily  all 
classes  of  flows  in  this  reference  before  it  can  be  accepted  as  a viable 
prediction  tool.  This  is  akin  to  a "calibration"  technique  for  TBLPM's. 

One  way  of  classifying  the  available  TBL  data  is  in  terms  of  the 
sign  of  the  applied  pressure  gradient  and  whether  or  not  the  flow  is  in 
equilibrium. 

Reference  [3]  hos  ranked  the  data  according  to  the  difficulty  en- 
countered by  the  28  calculation  methods  in  predicting  the  flows.  In 
order  of  increasing  difficulty,  these  were 

(a)  zero  and  mild  favorable  pressure  gradients, 

(b)  strong  favorable  and  mild  adverse  pressure  gradients, 

(c)  separating,  relaxing,  and  reattaching  flows. 

All  the  data  were  checked  for  two-dimensionality  by  normalizing  the 
momentum  integral  equation  and  integrating  in  the  x direction,  giving 


The  subscript  o indicates  quantities  at  the  start  of  the  flow. 

The  left  and  right  sides  of  this  equation  were  called  PL  and  PR,  respec- 
tively. If  the  measured  values  are  exact,  the  b.l.  two-dimensional,  and 
the  normal  stress  terms  negligible,  then  PL  = PR.  Since  all  these  con- 
ditions can  never  be  met  in  practice,  about  all  that  can  be  said  is  that 
PL  - PR,  and  that  strong  departures  from  this  equality  suggest  that 
some  or  all  of  these  conditions  are  not  met. 

Interestingly,  the  degree  of  difficulty  in  predicting  a flow  is 
directly  proportional  to  the  imbalance  between  PL  and  PR.  The  obvious 
conclusion  is  that  if  the  data  do  not  satisfy  the  2-D  momentum  integral 
equation,  then  a calculation  method  using  this  equation  cannot  predict 
the  data! 


C.  Three-Dimensional  Correction 

Assuming  that  the  imbalance  in  PL  and  PR  is  due  to  sidewall  b.l. 
growth  and  that  it  can  be  modeled  as  a source  or  sink  placed  along  a 
plane  of  symmetry,  Schlichting  [2]  showed  that  the  momentum  integral 
equation  can  be  balanced  by  including  a crossflow  term. 


^ + (2+H)  ^ -3— 
dx  U dx 


(3-2) 


where  x is  the  location  of  the  fictitious  source  or  sink,  and  may  be 
c 

obtained  by  solving  Eqn.  (3-2)  for  x^ , giving 


x 

c 


eu2/(0U2) 


x + 


- (PL-PR) 


(3-3) 


Unfortunately,  the  PL  and  PR  values  are  quite  noisy,  leading  to 
violent  fluctuations  in  the  value  of  xc . A second  method  is  to  arbi- 
trarily adjust  x^  to  give  the  best  match  with  experiments. 

There  are  serious  objections  to  using  this  correction  term.  Both 
methods  for  obtaining  x^  depend  on  having  experimental  data  available. 
This  can  hardly  qualify  as  a prediction  method!  This  correction  term 
will  not  be  used  for  a priori  diffuser  calculations.  However,  it  will 
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be  employed  in  this  chapter  during  the  calibration  process,  so  as  to 
enable  comparison  with  TBLPM's  in  [3], 

D.  Initial  and  Boundary  Conditions 

* 

Values  of  H,  6 are  known  at  the  starting  location,  and  the  ini- 
tial Tmax/,P  is  calculated  assuming  equilibrium  conditions.  The  pres- 
sure gradient  is  prescribed  and  6,  u^,  u^  calculated  along  the  flow. 
Instead  of  imposing  smoothed  dp/dx  values,  as  was  done  in  [3],  a ten- 
sioned spline  [48]  was  fitted  through  the  data  and  the  derivative  ob- 
tained numerically.  Three-dimensional  correct'  ms  were  applied  to 
detaching  flows  only. 

It  is  again  emphasized  that  detaching  diffuser  flows  have  to  be 
calculated  in  simultaneous  iteration,  and  the  present  prescribed  pres- 
sure gradient  mode  is  for  comparison  purposes  only. 

E.  Comparison  with  Experiments 

Predictions  using  the  UIM  equations  are  shown  in  Figs.  9 through  19. 
The  numbers  in  parentheses  after  each  flow  is  the  identification  assigned 
in  [3]. 

Close  agreement  is  obtained  with  data  for  accelerating  and  deceler- 
ating flows,  including  equilibrium  b.l.'s  of  both  types  (see  Figs.  9-13). 
Figs.  9a  and  9b  show  the  prediction  for  the  relaxing  flow  of  Bradshaw 

et  al.  (2400).  With  the  attached  lag  parameter  X = .025,  excellent 
* ** 

agreement  of  H,  6 , and  C^/2  are  obtained.  Fig.  9b  shows  the  devel- 
opment of  t /p  along  the  flow.  If  the  data  satisfied  Bradshaw's 
max 

correlation  exactly,  entrainment  and  maximum  shear  would  coincide  at 
every  station.  The  predicted  entrainment  rate  is  somewhat  low  at  the 
beginning,  but  is  quite  good  for  the  rest  of  the  flow.  The  low  starting 
value  is  probably  due  to  the  assumption  of  equilibrium  starting  condi- 
tions, while  in  fact  the  flow  is  far  removed  from  this  state. 

The  program  has  problems  predicting  the  reattaching  Tillmann  ledge 
flow  (1500),  Fig.  14.  H and  6*  are  overpredicted,  while  the  skin 
friction  values  are  too  low.  Nash  and  Hicks  [3]  were  able  to  improve 
agreement  with  data  by  doubling  the  initial  shear  stress  value.  This 
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is  expected  to  improve  the  present  prediction,  but  has  not  been  attempted 
since  this  type  of  flow  would  normally  be  calculated  in  simultaneous  it- 
eration. 


Problems  encountered  in  calculating  detaching  b.l.'s  with  prescribed 
pressure  gradient  (cyclic  iteration)  were  discussed  in  Chapter  One. 

These  are  evident  in  Figs.  15  through  18.  In  all  cases  the  flow  proceeds 
towards  detachment  but  relaxes  prematurely.  A 3-D  correction  term  was 
included  and  adjusted  to  give  the  best  possible  match  with  data.  The 
agreement  is  improved,  but  premature  detachment  occurs  if  too  small  a 
value  of  xc  is  used.  Moses'  asymmetric  diffuser  flow.  Fig  17,  and 
Strickland-Simpson's  airfoil  flow.  Fig.  18,  will  be  recalculated  with  a 
full  2-D  core  in  simultaneous  iteration  to  show  the  improved  prediction 
possible  (see  Chapter  Six). 

For  contrast.  Perry's  flow  (2900)  was  recalculated  as  a diffuser  with 
a 1-D  core,  as  described  in  the  next  chapter.  The  results.  Fig.  16,  bear 
out  the  claims  made  for  the  simultaneous  iteration  concept.  The  predicted 
b.l.  quantities  found  by  simultaneous  iteration  are  much  closer  to  the 
data  than  those  obtained  for  the  prescribed  pressure  gradient  method, 
especially  in  the  detaching  region.  In  calculating  the  flow,  it  was 
assumed  that  the  upper  and  lower  b.l.'s  were  identical  at  the  starting 
section.  This  is  not  the  case,  and  it  is  expected  that  improved  agree- 
ment would  result  if  the  correct  starting  values  were  available. 

Finally,  So  and  Mellor's  [46]  b.l.  growing  over  a convex  surface  is 
shown  in  Fig.  19.  The  results  are  in  accordance  with  Bradshaw's  [47] 
observation  that  the  turbulence  production  is  suppressed  on  a convex  wall 
b.l.  and  enhanced  on  a concave  one.  The  skin  friction  falls  along  the 
flow  as  the  turbulence  level  decreases,  and  the  prediction  is  too  high. 
This  flow  was  included  since  the  curved  throat  region  of  many  diffusers 
is  a convex  surface,  even  though  the  flow  length  is  small.  It  is  pres- 
ently not  known  how  long  the  curved  region  has  to  be  in  order  to  have  a 
significant  effect  on  the  downstream  flow,  nor  is  the  magnitude  of  the 
curvature  effect  well  established  at  this  time  (1976).  However,  this 
phenomenon  is  known  to  be  important  in  many  passage  applications,  and  the 
effect  needs  to  be  pointed  out  so  that  improved  TBLPM's  that  properly 
account  for  curvature  effects  will  be  created. 
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F.  Chapter  Conclusions 

(a)  The  current  TBLPM  is  capable  of  accurately  predicting  equilib- 
rium flows,  as  well  as  accelerating  and  decelerating  b.l.'s.  For  at- 
tached b.l.'s,  its  performance  is  as  good  as  the  best  prediction  method 
presented  in  the  1968  Stanford  Conference.  It  has  the  additional  advan- 
tage that  it  is  capable  of  calculating  detached  flows. 

(b)  The  method  has  no  difficulty  in  predicting  accelerating,  mildly 
decelerating,  and  equilibrium  flows.  For  detaching  flows,  the  inclusion 
of  the  3-D  correction  term  improves  the  accuracy  until  the  flow  nears 
detachment;  after  this  point  the  computed  values  are  no  longer  accurate. 

Inclusion  of  the  shear-stress  lag  equation  is  believed  to  be  the  reason 
for  the  good  prediction  of  strongly  perturbed  flows.  An  exception  is 
Tillmann's  reattaching  flow,  which  was  not  well  predicted.  B.l.'s  over 
curved  surfaces  are  not  well  predicted  either. 

(c)  The  procedure  does  extremely  well  for  b.l.'s  encountered  in  a 
typical  diffuser,  which  exhibit  mild  acceleration  in  the  inlet,  strong 
acceleration  around  the  throat,  and  strong  deceleration  in  the  diffusing 
section.  Thus  this  method,  when  used  in  cyclic  iteration  (prescribed 
pressure  gradient),  shows  the  weaknesses  seen  in  all  the  methods  of  the 
1968  Conference  for  flows  nearing  detachment.  However,  all  the  methods 
in  the  1968  Conference  also  use  cyclic  iteration.  As  shown  in  the  next 
chapter,  these  difficulties  near  detachment  do  not  occur  when  the  simul- 
taneous iteration  procedure  is  used. 
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CHAPTER  FOUR 

RESULTS  - ONE-DIMENSIONAL  CORE  DIFFUSERS 

A.  Discussion  of  Available  Data 

We  restrict  ourselves  to  two-dimensional  diffusers  and  briefly  re- 
view the  currently  available  data. 

Diffusers  have  either  straight  or  curved  centerlines,  as  depicted 
in  Fig.  K.  They  may  be  symmetric  or  asymmetric  about  the  centerline. 

The  most  widely  studied  is  type  (a),  for  which  flow-regime  charts 
were  established  by  Fox  and  Kline  [1].  Reneau  et  al.  [45]  created  a set 
of  data  maps  that  can  be  used  to  estimate  the  overall  pressure  recovery, 
Cp.  Carlson  et  al.  [27]  compared  the  performance  of  types  (a)  and  (b). 

Fox  and  Kline[30]  established  flow-regime  charts  for  type  (c),  sometimes 
called  a circular  arc  diffuser.  Sagi  et  al.  [31]  made  measurements  of 
both  types  (c)  and  (d). 

In  all  the  above  experiments,  only  "zeroth-order"  quantities  were 
measured.  These  were  Cp(x)  and  the  gross  qualitative  features  of  the 
flow,  such  as  the  levels  of  unsteadiness  from  visualization  of  wall 
tufts  and  whether  or  not  backflow  was  present  in  an  intermittent  or 
steady  basis.  The  b.l.  integral  thicknesses  at  the  inlet  were  recorded. 
There  were  no  detailed  measurements  of  the  b.l.  development  along  the 
flow,  and  no  skin  friction  or  turbulence  data  were  taken. 

Moses  [7]  measured  the  variation  of  Cp(x)  and  integral  parameters 
along  the  wall  of  a type  (e)  diffuser  in  transitory  stall.  Unfortunately, 
the  diffuser  throat  had  a small  radius,  and  it  is  possible  that  strong 
curvature  effects  may  have  introduced  unexpected  behavior  in  the  b.l. 

The  most  extensive  data  for  a single  unit  available  today  is  the 
airfoil  type  flow  of  Strickland  and  Simpson  [32],  also  on  a type  (e) 
diffuser.  Detailed  measurements  of  the  b.l.  development  are  available, 
including  details  of  the  turbulence  quantities  along  the  flow.  These 
data  were  taken  with  a directionally  sensitive  laser  anemometer,  so  that 
the  measurements  are  expected  to  be  more  accurate  than  pitot  or  hot-wire 
data  in  regions  where  the  fluctuations  are  large  and  the  meanflow  direc- 
tion uncertain. 
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B . Results  and  Discussion 

Figures  20  through  27  compare  predictions  of  the  current  method  with 

data.  Figs.  20  and  21  compare  the  variation  in  velocity  ratio  u/u__„ 

REF 

along  the  diffuser  walls  with  the  data  of  Carlson  et  al.  [27].  Fig.  20 
is  for  an  unstalled  diffuser,  for  which  the  velocity  variation  is  pre- 
dicted within  the  uncertainty  of  the  data.  Fig.  21  is  for  a diffuser  in 
transitory  stall,  and  the  calculated  values  are  barely  outside  the  uncer- 
tainty band  in  the  region  between  the  intermittent  and  fully  developed 

* 

detachment  points.  Fig.  22  shows  the  predicted  Variation  in  H,  6 , 
and  C^./2  along  this  diffuser  — no  data  are  available  for  comparison. 

The  predictions  for  a complete  series  of  N/Wl  = 6,  B1  = .030  dif- 
fusers with  area  ratios  (AR)  varying  from  1.5  to  2.65  are  shown  in  Fig. 
23.  The  calculations  compare  very  well  with  the  measurements  of  Carlson 
et  al.  [27].  The  Cp  values  from  the  data  maps  of  Reneau  [45]  are  some- 
what higher,  but  are  well  within  the  uncertainty  of  the  data. 

The  predicted  exit  conditions  for  the  same  series  of  diffusers  are 
shown  in  Fig.  24.  Only  one  data  point  is  available,  for  AR  = 1.8.  The 
agreement  in  this  case  is  excellent. 

Figure  25  is  a replot  of  the  predicted  variation  in  exit  Cp,  in 

addition  to  which  are  shown  the  locations  of  the  intermittent  (H  = !!„„„) 

SEP 

and  fully  developed  (C^/2  = 0)  detachment  points  as  fractions  of  the 
diffuser  length  (x/L) . Also  shown  are  the  locations  designated  TI 
(incipient  transitory  stall)  and  IT  (intermittent  transitory  stall)  from 
flow  visualization  of  Cajrlson  et  al.  [27].  For  the  few  data  points 
available,  the  TI  location  is  quite  well  predicted  by  the  intermittent 
detachment  point  according  to  the  Sandborn  criterion.  The  location  of 
fully  developed  detachment  occurs  somewhat  downstream  of  this  point. 

Figure  26  is  a summary  of  the  performance  of  N/Wl  = 12  diffusers 
as  a function  of  the  divergence  angle  29,  for  the  inlet  blockage  B1 
varying  from  .007  to  .050.  The  accuracy  of  the  prediction  decreases  as 
20  and  B1  increase.  This  conclusion  is  in  accordance  with  the  find- 
ings of  Woolley  et  al.  [4],  For  small  Bl,  the  b.l.  is  a correction  to 

★ 

the  throughflow,  so  that  small  errors  in  6 cause  even  smaller  errors 
in  Cp.  As  the  b.l.  becomes  a significant  portion  of  the  flow,  the  accu- 
racy of  the  Cp  predictions  decreases.  The  data  for  Bl  = .050  have  a 
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sharp  peak,  following  which  it  levels  off  at  a constant  value  of  Cp. 


The  beginnings  of  this  peaky  behavior  can  be  seen  in  the  curve  for  B1  = 

.030.  The  calculation  is  unable  to  follow  this  trend.  The  behavior  of 
the  calculated  results  is  similar  for  all  inlet  blockages. 

The  locations  marked  X indicate  the  value  of  20  at  which  the 
shear  layers  from  the  upper  and  lower  walls  begin  to  interfere  with  each 
other.  No  irrotational  region  remains,  and,  viewed  strictly,  the  calcu- 
lation method  is  not  valid  beyond  this  point. 

Finally,  Fig.  27  is  a summary  of  all  diffusers  that  have  been  run. 

The  current  method  is  able  to  predict  Cp  of  all  tested  units  to  about 
the  uncertainty  in  the  data,  with  the  exception  of  the  B1  = .050  case. 

For  all  Bl,  the  present  calculation  is  capable  of  predicting  Cp 
within  ± 6%  of  the  data  for  all  diffusers  whose  divergence  angle  29  is 
less  than  1.2  x 20  . The  range  of  calculation  has  therefore  been 

SL~  SL 

doubled  from  20/28  =0.6  in  the  method  of  Woolley  and  Kline  [4]  to 

29/20  =1.2  in  the  present  method.  This  extension  carries  the  method 

a a * 

well  into  a region  beyond  the  peak  Cp  — up  to  approximately  the  line 
of  appreciable  stall. 

C.  Chapter  Conclusions 

(a)  Tue  diffuser  calculation  method  assuming  a 1-D  core  and  symmet- 
rical b.l.'s  is  capable  of  predicting  the  performance  in  the  transitory 
stall  regime  well  past  the  peak  in  the  Cp  curve.  Accuracies  of  ± 6% 
can  be  obtained  up  to  divergence  angles  that  are  1.2  times  the  location 

of  20  , even  when  the  difficult  case  of  Bl  = .050  is  included, 

a-a 

(b)  Prediction  accuracy  decreases  with  increasing  inlet  blockage. 

Neglecting  the  highest  blockage  value  of  .050,  the  ± 6%  accuracy  in  Cp 

can  be  obtained  for  20/20  = 1.8. 

a-a 

(c)  The  predicted  location  of  intermittent  detachment  (H  = H ) 

b hr 

according  to  the  Sandborn  criteria  occurs  very  close  to  the  point  desig- 
nated as  "incipient  transitory  stall" (TI)  in  the  flow  visualization  data 
of  diffusers.  The  location  of  zero  wall  shear  occurs  a small  distance 
downstream  of  this  point. 

(d)  The  program  was  tailored  to  model  detached  flows  using  very 
sparse  data.  The  program  output  consists  of  turbulent  quantities,  wall 
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shear  stresses  and  entrainment  values  for  which  experimental  data  are 
almost  nonexistent.  Much  more  data  of  detailed  b.l.  development  and 
turbulence  quantities  in  the  detaching  regions  is  needed  to  extend  the 
applicability  and  either  fully  verify  or  improve  the  present  model. 
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CHAPTER  FIVE 

DEVELOPMENT  OF  PREDICTION  METHOD  FOR  2-D  CORE  DIFFUSERS 

A.  Limitations  of  the  1-D  Core  Method 

The  predictions  of  symmetric  diffusers  with  symmetric  b.l.'s  using 
the  1-D  irrotational  core  model  were  shown  to  be  quite  acceptable  for 
engineering  purposes. 

There  are  many  cases,  however,  for  which  such  a 1-D  core  is  ob- 
viously a poor  approximation,  such  as  for  a grossly  asymmetrical  diffu- 
ser. Not  so  obvious  is  the  fact  that  the  use  of  the  1-D  core  approxi- 
mation and  the  simultaneous  streamwise  marching  procedure  has  enabled  us 
to  convert  from  an  elliptic  problem  to  a fully  parabolic  one.  Certain 
essential  information  has  been  inevitably  lost  in  this  process.  The 
flow  of  information  in  the  numerical  procedure  is  in  the  downstream  di- 
rection only  — all  upstream  propagation  is  totally  absent. 

It  is  well  known  that  the  effect  of  downstream  blockage  can  play  an 
important  role  in  determining  the  upstream  pressure  gradient  and  hence 
can  control  the  detachment  process.  This  elliptic  field  effect  can  be 
clearly  seen  in  the  experiment  of  Chui  et  al.  [33]  on  a fully  stalled 
diffuser.  The  dominant  adverse  pressure  gradient  occurs  well  ahead  of 
the  diffuser  throat,  and  is  due  mainly  to  the  blockage  of  the  stalled 
flow  downstream  of  the  throat.  The  elliptic  field  effect  causes  the 
streamlines  to  diverge  ahead  of  the  throat,  in  a region  that  is  bounded 
by  parallel  walls.  A 1-D  calculation  would  predict  acceleration  in  this 
region  and  could  obviously  never  predict  this  flow  even  approximately 
for  the  lowest-order  quantities. 

The  importance  of  the  elliptic  field  effect  in  the  transitory  stall 
regime  is  not  known.  It  is  probably  not  as  pronounced  as  in  the  fully 
stalled  regime,  on  account  of  the  smaller  detachment  zones  and  the  con- 
sequent smaller  curvatures  of  the  streamlines. 

There  is  a basic  dilemma,  however.  It  was  pointed  out  in  Chapter 
One  that  detached  b.l.'s  can  only  be  calculated  simultaneously  with  the 
adjacent  Irrotational  freestream.  Numerical  stability  requires  that  the 
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pressure  gradient  acting  on  a b.l.  at  and  near  detachment  be  exactly  that 
which  occurs  as  the  result  of  mutual  interaction  between  it  and  the  irro- 
tational  core.  The  solution  to  Laplace's  equation  in  the  core  requires 
information  along  the  entire  boundary  in  order  to  be  well  posed.  There 
is  thus  a basic  conflict  between  the  requirements  of  the  b.l.  and  the 
inviscid  core.  The  final  solution  will  therefore  have  to  be  obtained 
iteratively,  each  iteration  being  designed  in  such  a manner  as  to  satisfy 
these  separate  and  conflicting  requirements. 

B . Procedure  for  Calculation  of  Diffusers  with  2-D  Core 

The  basic  outline  of  the  calculation  method  for  diffusers  with  2-D 
core  will  now  be  developed.  The  next  few  sections  present  the  b.l.  and 
potential  flow  schemes. 

In  cyclic  iterative  calculations  of  the  type  used  by  Woolley  [4]  and 
White  [5],  the  solution  of  Laplace's  equation  in  the  domain  bounded  by 
the  diffuser  walls  gives  an  estimate  of  the  velocity  gradient  in  the 
streamwise  direction,  which  is  prescribed  to  calculate  the  b.l.  growth. 

k 

The  displacement  thickness  6 is  subtracted  from  the  diffuser  walls  to 
give  a new  effective  flow  channel  (EFC) . Laplace's  equation  is  solved 
in  this  new  EFC,  giving  a new  estimate  of  the  velocity  gradient,  which 

•k 

is  used  to  obtain  a new  6 , etc.  The  process  is  considered  to  have 

* 

converged  if  the  change  in  6 or  the  velocity  gradient  between  success- 
ive iterations  is  smaller  than  a preselected  tolerance. 

This  scheme  works  for  unstalled  diffusers  and  for  the  fully  stalled 
case  for  which  the  simplifying  assumption  of  constant  pressure  in  the 
stalled  zone  permits  the  detached  b.l.  to  be  modeled  as  a free  stream- 
line problem.  In  this  case  of  a fully  stalled  diffuser,  the  prescribed 
pressure  gradient  calculation  is  terminated  before  intermittent  detach- 

k 

ment,  and  the  6 line  is  extrapolated  into  the  stalled  zone,  where  its 
location  is  iteratively  determined.  For  detaching  b.l.'s,  however,  cyc- 
lic iteration  is  numerically  unstable,  and  this  procedure  diverges. 
Further,  the  pressure  in  the  transitory  stalled  zone  is  not  constant, 
and  substantial  pressure  recovery  occurs  in  this  state,  so  that  a free 
streamline  model  is  inappropriate.  A new  scheme  that  avoids  these  prob- 
lems is  needed. 
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An  examination  of  the  data  of  Moses  [7]  and  Strickland  et  al.  [32] 
shows  that  the  velocity  profiles  in  the  irrotational  core  of  diffusers 
operating  in  the  transitory  stall  regime  are  quite  linear,  excluding,  of 

: k 

course,  regions  of  sharp  curvature  in  6 such  as  near  the  throat. 

That  is,  data  show  that  there  is  a linear  variation  in  edge  velocity  be- 
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tween  the  upper  and  lower  6 lines.  In  fact,  after  detachment,  the 
profile  becomes  almost  one-dimensional,  which  may  be  why  the  1-D  core 
method  is  so  successful. 

We  note  further  from  the  data  of  Smith  and  Kline  [34]  that  transi- 
tory stall  begins  and  is  restricted  to  one  wall,  even  if  the  two  walls 
of  the  diffuser  are  nominally  symmetric.  This  is  not  surprising  since 
when  one  b.l.  detaches,  the  pressure  gradient  on  the  opposite  wall  is 
immediately  relaxed.  In  an  asymmetric  diffuser  there  is  no  question 
that  detachment  is  restricted  to  the  diverging  wall. 

The  diffuser  will  therefore  be  modeled  with  an  upp-  r wall  that  has 
an  attached  b.l.  and  a lower  wall  in  which  the  layer  may  or  may  not  be 
detached.  The  upper  b.l.  and  that  portion  of  the  lower  b.l.  well  ahead 
of  detachment  can  be  calculated  with  prescribed  pressure  gradient,  while 
the  detaching  and  detached  regions  have  to  be  simultaneously  calculated 
with  the  inviscid  core. 

If  the  edge  velocity,  from  the  lower  wall  b.l.  calculation 

is  the  same  as  the  edge  velocity,  obtained  from  the  solution  of 

Laplace's  equation  in  the  EFC  defined  by  the  new  6 lines,  the  solu- 

tion is  considered  to  have  converged.  Otherwise,  the  process  is  contin- 

* 

ued  by  prescribing  on  the  upper  wall  and  using  the  6^  so  ob- 

tained to  define  a new  EFC.  Laplace's  equation  solved  in  this  new 
domain  gives  a new  against  which  the  edge  velocity  from  a 

new  lower  wall  b.l.  calculation  may  be  compared,  and  so  on. 

The  convergence  criterion  used  is  that  for  all  stations,  Cp  <_ 
Cperor,  where 


Cperor 

= 

|CplD  " Cp2D | 

where 

CplD 

1 " ^U'°1D^UREF 

and 

Cp2D 

= 1 - (U»2D/UREF 
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A Cperor  = .025  can  be  achieved  in  8 to  10  iterations,  and  has 
been  used  for  the  predictions  described  in  the  next  chapter. 

The  only  regions  where  the  linear  variation  in  velocity  between  the 
upper  and  lower  walls  breaks  down,  is  in  those  areas  where  the  streamlines 
are  sharply  curved,  such  as  near  the  throat  and  in  the  rapidly  growing 
zone  ahead  of  intermittent  detachment.  The  b.l.  is  accelerating  around 
the  throat  and  can  therefore  be  calculated  with  prescribed  pressure  gra- 
dient. The  region  of  strong  streamline  curvature  is  also  well  ahead  of 
detachment,  and  it,  too,  can  be  calculated  in  a similar  manner.  The 
lower  wall  is  therefore  calculated  with  prescribed  pressure  gradient  and 
switched  over  to  the  simultaneous  linear  velocity  profile  scheme  for 
H > 0.9  where  H is  the  Sandborn  criterion.  The  entire  pjroc- 

— bcr  obr 

ess  is  shown  in  Figs.  L and  M. 

In  summary,  the  present  scheme  is  broadly  similar  to  a predictor- 

corrector  method.  The  linear  core  profile  method  is  the  predictor,  which 

* 

provides  an  estimate  of  the  lower  wall  6 and  edge  velocity  U , „ . 

s & 3 ^lDs 

The  corrector  is  the  values  of  the  edge  velocity  U obtained  from 

the  solution  of  Laplace's  equation  in  the  EFC.  When  the  predicted  and 
corrected  values  agree  within  an  acceptable  tolerance,  a converged 

solution  is  obtained.  The  final  solution  reflects  the  accuracy  of  the 
corrector,  and  the  approximations  of  the  predictor  are  no  longer  present. 
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A = prescribed  pressure  gradient  calculation 
B = simultaneous  linear  core  profile  method 


Fig.  L. 


Sketch  of  the  2-D  core  diffuser  illustrating  regions 
where  the  two  types  of  calculation  methods  are  used. 


C . Simultaneous  B.L.  Calculation  with  Linear  Core  Velocity  Profile 

k 

The  location  of  the  upper  wall  <$^  line  is  known  from  the  last 

iteration,  as  is  the  velocity  distribution  U from  the  2-D  potential 

°°2Du 

flow  calculation  in  the  resulting  EFC.  The  diffuser  width  W is  known 

k 

from  the  input  geometry.  We  wish  to  calculate  the  lower  wall  and 

the  corresponding  edge  velocity,  U , assuming  linear  variation  of 
velocity  between  the  upper  and  lower  6 lines.  The  figure  below  shows 
the  situation  at  location  x. 
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START 


Solve  V <J>  = 0 with  bare 
diffuser  walls  as  EFC. 


Estimate  6u  from  flat  plate  solution. 


Assume  linear  core  profile  aijd  calculate  lower  wall 

b.l.  in  channel  defined  by  5 and  bottom  wall. 

u 


Locate  new  EFC  bounded  by 

6*  and  6*.  Solve  v2(j>  ■ 0. 
u s 


Is  |Cp.n  - CP-_ | < Cp^ 


1 *1D  r2D 1 error  \ Yes 
in  entire  channel?  / 


Is  this  an  \ No 
even  iteration?  / 


Prescribe  U—_n  011 

upper  wall  b.l.  and 

calculate  6 . 

u 


' Is  \ Yes 
vX  - XEXIT?/ 


Continue  lower  b.l.  calculation 
using  linear  core  profile  in 
channel  defined  by  6*  and  ^ 
bare  bottom  wall.  Calculate  6 . 


• 

^ i ^ 

Fig.  M.  Flowchart  illustrating  the  two-dimensional  core  calculation 
method . 
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A 

Knowns:  U _,,  , 6 , W 

“2Du  u 

* 

Unknowns:  U , 6 
“IDs  s 


The  volumetric  flow  at  section 


■ (“  - < - 5‘)p 


“2Du+^“lDs  ^ 


(5-2) 


Define 


* * 

6-6 
u s 


U on  + U in 

“2Du  “IDs 


Differentiating  Eqn.  (5-2)  with  respect  to 
from  mass  conservation,  and  rearranging  gives 


setting  dQ/dx  = 0 


(=$1 


m\4^\ 

oo  oo 

1=  (w  - <)! 


In  the  above  equation,  U 


has  been  written  as  U for  brevity. 


The  right-hand  side  of  this  equation  is  known  from  previous  iteration. 

Therefore,  Eqn.  (5-5)  can  be  used  to  replace  the  1-D  core  equation  (2-43) 

and  the  new  set  of  equations,  (2-26),  (2-32),  (2-34),  and  (5-5)  solved 

in  a stepwise  fashion  along  the  flow. 

The  dependent  variables  can  be  processed  as  before  to  obtain  the 
•k  k 

location  of  the  lower  6g  line,  which,  together  with  the  upper  6^  line 
obtained  in  the  last  iteration,  forms  the  boundary  of  the  EFC. 

The  2-D  Laplace  solver  is  next  outlined. 

D.  Solution  of  the  2-D  Laplace  Equation 

We  desire  to  solve  Laplace's  equation  in  the  domain  bounded  by  the 

★ 

upper  and  lower  6 lines  and  the  inlet  and  exit  planes  of  the  diffuser. 
The  velocity  is  assumed  constant  across  the  inlet,  and  it  is  specified 

k 

that  there  is  no  flow  across  the  upper  and  lower  6 lines,  which  are 
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thus  approximated  as  streamlines.  The  situation  is  depicted  in  the 
following  figure. 


ufl02DU 


V2<*>  = 0 IN  EFC 


U©2DS 


JBS* 


The  edge  velocity  U „ is  needed  along  the  entire  boundary  of  the 
EFC.  This  is  similar  to  the  problem  solved  by  Woolley  et  al.  [4],  and 
lends  itself  naturally  to  a boundary  integral  method,  since  only  the 
values  along  the  boundary  are  required. 

Two  shortcomings  of  the  method  used  in  [4]  were  that  the  exit  veloc- 
ity was  assumed  to  be  one-dimensional  and  the  equation  formulation  used 
the  Cauchy-Riemann  conditions,  which  necessitated  the  taking  of  numerical 
derivatives  with  their  potential  for  large  errors. 

Recently  my  colleague  Rinehart  [36]  has  developed  a similar  method 
for  solving  the  2-D  Laplace  equation  which  avoids  both  these  difficul- 
ties. Since  his  work  is  soon  to  be  published,  only  an  outline  of  the 
method  will  be  presented. 

Consider  a simply  connected  domain  V in  the  complex  plane  bounded 
by  a smooth  closed  contour  C.  Let  z^  be  an  interior  point,  and,  if 
f(z)  is  analytic  in  V,  Cauchy's  integral  formula  gives  the  value  of 
the  function  at  this  point  as 


27Ti  f ( z ) 
o 


f z-z 
c o 


(5-6) 


Now  let  z approach  the  contour  C.  In  the  limit  when  z is 
o o 

on  C,  we  have  the  Plemelj  formula. 


Itt  f(z  ) ■ P f dz 

o J z-z 

n r\ 


(5-7) 


■ 


I 

f 


: ’-T»  jr 


KJLil 


■ 


j 


k 

I 


The  integral  on  the  right-hand  side  is  to  be  interpreted  in  the 
Cauchy  principal  value  sense. 


If  C is  not  a smooth  curve  and  zq  is  a corner  point,  Eqn.  (5-7) 
is  modified  to 


ia  f(z  ) = P / dz  (5_8) 

o J z-z 

c o 

where  a is  the  interior  angle  at  the  corner.  For  a smooth  curve, 
a = tt  and  Eqn.  (5-8)  reduces  to  (5-7).  For  further  details,  see  Muskhe- 
1 ishvi 1 i [ 37 ] . 

The  boundary  C is  discretized  into  N segments  whose  end  points 
are  numbered  increasing  in  the  counterclockwise  direction,  as  shown  on 
the  next  page. 

Let  the  boundary  point  zQ  at  which  the  function  is  to  be  evalua- 
ted be  located  at  node  C . Then,  since  the  singularity  is  present  at 

m 

"this  point  alone,  Eqn.  (5-7)  can  be  rewritten  as  the  sum  of  an  ordinary 
contour  integral  plus  a principal  value  integral. 


iff  f ( z ) 
o 


p f 

J °PV  Z~Zo 


f (z) 

z-z 

o 


dz 


(5-9) 


f(z)  is  expanded  in  a separate  Taylor  series  expansion  along  each  in- 
terval of  the  boundary,  and  on  performing  integrations  of  the  resulting 
terms  we  get 
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1 


m-2 


m+1 


in  f (z 


0 k>  = £ M + £ f,A  , 

°’k  j=m+l  J Jk  j=m-l  J JK 


k = 1,  2,  ...  , M-l 

(5-10) 


where 


f (Zj) » 


Jk 


the  geometry  factors  for  the  principal  value  segment  at  the 
node  j when  the  singularity  is  at  node  k. 


Jk 


= the  geometry  factors  of  the  rest  of  the  boundary. 


Cpn-1  ^m+l 


CPV,  PRINCIPAL 
VALUE  SEGMENT 


th 


Rewriting  Eqn. 

(5-10)  for  the  m node 

m+N-2 

P P 

j&i  f)Ajk  + 

f , , , + f (Ar. -Itt)  + 

m-l  m- 1 , k m mk 

II 

u 

o 

Def ine 

Ajk  * Gjk  + iHjk 

P p p 

Ajk  = Gjk  + iHjk 

= 0 , 


m+1  m+l,k 
for  k ■ 1,  2,  ...  , (N-l)  . (5-11) 


(5-12) 
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Choosing  the  analytic  function  as  f(z)  = 2n  V - ia,  where  V is 
the  magnitude  of  the  velocity  and  a is  the  local  streamline  angle. 


a = tan 


" (*)• 


we  have 


fj  = *n  v - io  , j = 1,  2 N . (5-13) 

Substituting  Eqns.  (5-12)  and  (5-13)  into  (5-11)  and  taking  the 
imaginary  part  gives 


nri-N-2 


m+N-2 


E *n  vi  Im(Aik)  + D *n  vi  Zm(Aik)  - n £n  vk  = E «4*, (a ik> 

j =m+l  J JK  j =in- 1 J JK  * j=m+l  J 6 


urrx 

+ V ct.R  (A ) , 

. *-> , j e'  jk'  ’ 

j=m-l  J J 


k - 1,  2,  ...  , (N-l) 


(5-14) 


This  set  of  equations  may  be  written  as  the  matrix  equation, 

A i.n  V = b,  as  shown  on  the  next  page. 

p 

The  geometry  factors  Aj^  and  A ^ ^ are  all  known,  as  are  the  . 
Eqn.  (5-15)  can  thus  be  solved  for  the  unknown  , the  velocities  along 
the  EFC. 

Given  the  geometry  of  the  EFC,  the  velocities  along  its  edge  can 
thus  be  calculated. 
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CHAPTER  SIX 

RESULTS  - TWO-DIMENSIONAL  CORE  DIFFUSERS 

A.  Moses1  Asymmetric  Diffuser  Flow 

Moses'  diffuser  was  of  type  (e)  with  one  diverging  wall  at  an  angle 
of  11.31  degrees,  AR  = 2.5,  L/Wl  = 7.5,  and  the  b.l.  thickness  at  the 
throat,  0/W1  = .007.  The  sharp  throat  radius,  R/Wl  = .57,  caused  con- 
vergence problems  because  of  the  rapid  change  of  Cp(x)  in  the  throat 
region.  An  artificial  increase  of  R/Wl  to  1.0  allowed  convergence  with- 
out materially  affecting  the  downstream  solution.  A 3-D  correction  with 
Xc  = 100.0  ft  was  necessary  to  match  the  data.  The  results  are  shown  in 
Figs.  28  and  29. 

Cp  on  both  walls  is  predicted  to  the  accuracy  in  the  data,  which  is 

estimated  to  be  ± 6%.  The  qualitative  trends  of  Cp(x)  are  correct, 

including  the  sharp  suction  peak  at  the  throat  of  the  diverging  wall, 

and  the  steady  increase  on  the  unstalled  wall.  The  suction  peak  value 

is  considerably  underpredicted,  but  the  data  here  are  quite  questionable 

on  account  of  the  rapid  streamwise  variation  of  Cp  in  this  region. 

The  greatest  deviation  from  the  data  occurs  in  the  region  of  detachment. 

* 

H and  6 are  quite  well  predicted  before  detachment,  but  are  consid- 
erably overpredicted  in  the  reversed  flew  region. 

B.  Strlckland-Simpson  Airfoil  Type  Flow 

This  flow  is  in  a type  (e)  diffuser  with  a flat  bottom  wall.  The 
top  wall  converges  and  then  diverges,  giving  a pressure  distribution 
similar  to  that  on  the  upper  surface  of  an  airfoil. 

The  flow  was  calculated  with  prescribed  pressure  gradient  up  to 
8.11  ft,  at  which  point  the  experimenters  had  to  remove  most  of  the  upper 
wall  b.l.  to  force  the  flow  to  detach  on  the  lower  wall.  The  rest  of  the 
flow  was  calculated  simultaneously  with  a full  2-D  inviscid  core. 

Initial  attempts  to  predict  this  flow  resulted  in  very  rapid  growth 

* 

of  <5  and  H after  detachment,  similar  to  that  for  the  Moses  diffuser 
flow.  To  prevent  this  through  a 3-D  correction  would  have  required 
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negative  values  of  , which  is  not  realistic  for  a decelerating  flow 

with  growing  sidewall  b.l.'s.  Instead,  the  lag  equation  was  removed 

after  detachment,  giving  the  results  shown  in  Figs.  30  and  31. 

The  b.l.  growth,  H,  6 , and  C^/2  for  both  the  upper  and  lower 

walls  are  very  well  predicted,  except  for  a small  deviation  near  the 

exit.  The  location  of  both  intermittent  and  fully  developed  detachment 

is  closely  predicted,  but  the  skin-friction  values  in  the  reversed  flow 

region  are  somewhat  smaller  in  magnitude  than  the  data.  C^/2  for  the 

upper  wall  is  slightly  overpredicted,  but  the  uncertainties  in  these 

data  are  quite  large  on  account  of  the  thinness  of  this  b.l.  and  the 

consequent  poorer  definition  the  wall  regions  of  the  velocity  profiles. 

Figure  31  shows  the  variation  of  Cp(x),  U , and  the  nondimensional 

entrainment  rate  yr  dQ/dx  . U is  underpredicted  by  about  5%  in  the 

detached  region,  leading  to  a 6%  overprediction  of  Cp.  The  entrainment 

rate  is  quite  good  until  detachment,  when  it  abruptly  rises  in  response 

to  the  removal  of  the  lag  equation.  The  value  is  almost  100%  too  large 

at  detachment,  following  which  the  deviation  begins  to  decrease.  The 

reason  for  the  excellent  agreement  of  the  mean  flow  parameters  using  this 

incorrect  value  of  the  entrainment  rate  is  not  known.  It  is  a peculiar 

2 

coincidence,  however,  that  the  values  of  t /pU  computed  from  the 

max  00 

data  using  Eqn.  (2-40)  and  plotted  in  Fig.  8 display  this  same  trend. 

The  maximum  shear  stress  computed  from  the  data  also  have  their  largest 
deviation  near  the  detachment  point. 

C.  Discussion 

Both  diffusers  used  for  comparing  with  the  2-D  core  calculation  are 
type  (e) , with  one  diverging  wall,  these  being  the  only  data  available. 
This  is  an  unfortunate  choice,  since  the  flow  regimes  for  asymmetric  dif- 
fusers are  expected  to  be  somewhat  different  from  those  for  symmetric 
units.  Since  the  divergence  is  limited  to  one  wall,  the  b.l.  on  this 
wall  begins  to  detach  much  earlier  than  on  a symmetric  unit  with  the  same 
28.  Line  a-a  therefore  occurs  at  a lower  28  and  the  entire  flow  re- 
gime shifts  downward. 

Preferential  stall  occurs  and  is  restricted  to  the  diverging  wall. 
The  transitory  stall  regime  is  expected  to  be  almost  nonexistent  for 
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asymmetric  diffusers,  the  flow  changing  from  an  essentially  unstalled  to 
a quasi-steady  fully  stalled  flow  as  the  divergence  angle  is  increased 
at  constant  L/Wl.  The  limited  data  available  support  this  description. 

As  a consequence,  both  diffusers  are  actually  operating  in  the 
fully  stalled  mode  with  a relatively  steady  recirculating  separation 
bubble,  even  though  they  should  both  be  in  the  small  transitory  stall 
regime,  according  to  the  flow  regime  chart,  Fig.  1.  The  present  calcu- 
lation method  was  not  designed  for,  and  does  not  give  accurate  values  for, 
b.l.  parameters  in  the  fully  stalled  zone,  even  though  the  zeroth-order 
quantities,  the  Cp,  and  locations  of  detachment  are  quite  well  predic- 
ted. The  justification  for  removal  of  the  lag  equation  in  the  reversed 
flow  region  is  that  the  detached  lag  parameter  A^  was  determined  by 
matching  data  from  a diffuser  operating  in  the  transitory  stall  regime, 
while  the  Strlckland-Simpson  flow  is  closer  to  that  of  a fully  stalled 
case.  It  appears  from  the  good  predictions  obtained  with  no  lag  equa- 
tion in  detached  flow,  that  perhaps  a higher  A^  is  appropriate  in  this 

zone,  since  A,  °°  corresponds  to  an  instantaneous  response  between 
d 

the  local  velocity  profile  and  the  shear  stresses. 

An  interesting  feature  of  the  Str ickland-Simpson  flow  is  the  region 
in  the  neighborhood  of  partial  removal  of  the  upper  b.l.  at  the  entrance 
to  the  diffusing  section.  The  upper  b.l.  undergoes  a severe  perturba- 
tion and  slowly  relaxes. 

The  largest  deviation  from  the  data  in  all  the  diffusers  that  have 
been  run  occurs  in  the  region  between  intermittent  detachment  and  the 
location  of  zero  wall  shear.  This  is  evident  in  Fig.  29  (2-D  Moses  dif- 
fuser) and  Figs.  20  and  21  (1-D  Carlson  diffuser).  The  present  calcula- 
tion evidently  cannot  model  the  flow  closely  in  this  region.  The  agree- 
ment improves  both  upstream  and  downstream  of  this  zone. 

The  reasons  for  this  deviation  may  be: 

(a)  The  Coles'  profile  does  not  adequately  represent  measured 
velocity  profiles  in  the  neighborhood  of  zero  wall  shear. 

(b)  The  eddy-viscosity  formulation,  Fig.  8,  has  the  greatest  devi- 
ation from  data  in  this  region. 

(c)  The  effect  of  neglected  terms  in  the  momentum  integral  equa- 
tion, such  as  the  normal  stress  terms,  is  greatest  in  the  detachment  zone. 

56 


mm 


1 

I 

(d)  The  turbulence  measurements  have  the  greatest  uncertainty  in 
this  region  on  account  of  the  small  mean  and  large  fluctuation  magnitudes. 
Considering  all  these  negative  factors,  the  overall  success  of  the 
» current  method  is  gratifying. 
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CHAPTER  SEVEN 
SUMMARY 


A. 


Conclusions 


(a)  A calculation  method  has  been  developed  that  successfully  pre- 
dicts three  types  of  flows: 

•••  turbulent  boundary  layers  with  prescribed  pressure  gradient, 

•••  symmetric  diffusers  with  one-dimensional  core, 

•••  diffusers  with  two-dimensional  inviscid  core. 

The  last  two  types  can  have  attached  or  detached  boundary  layers. 

(b)  Diffuser  predictions  to  about  ± 6%  accuracy  in  Cp  can  be 

made  up  to  about  the  location  of  the  line  of  appreciable  stall  in  the 
transitory  stall  regime.  This  corresponds  to  20/20^  = 1.2.  Predic- 


a-a 

tion  accuracy  increases  with  decreasing  inlet  blockage. 

•k 


(c)  The  mean  boundary  layer  parameters  H,  6 , C^/2,  etc.,  are 


extremely  well  predicted.  For  diffusers,  the  locations  of  both  inter- 
mittent detachment  and  zero  wall  shear  are  also  predicted  with  remarkable 
accuracy.  However,  skin  friction  and  entrainment  in  the  reversed  flow 
region  are  only  fair. 

(d)  Execution  times  for  the  program  on  an  IBM  370/168  are  on  the 
order  of  0.25  seconds  for  a straight  boundary  layer  calculation,  and  1.0 
sec  for  a 2-D  Laplace  equation  solution.  A 1-D  core  diffuser  takes  about 
0.5  sec.  A typical  full  2-D  calculation  involves  6 to  10  iterations  of 
the  boundary  layer  and  inviscid  core  and  takes  about  10  secs  to  execute. 

(e)  The  overall  success  of  the  method  legitimizes  the  concept  of 
simultaneous  iteration  as  a means  of  preventing  the  singular  behavior  of 
the  classical  boundary  layer  calculations  in  the  neighborhood  of  detach- 
ment . 

(f)  The  eddy-viscosity  scheme  used  in  this  report  was  based  on  ex- 
tremely sparse  information.  Improved  predictions  will  be  possible  only 
when  more  data  on  detached  and  detaching  boundary  layer  behavior  become 
available . 
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B.  Recommendations  for  Further  Work 


(a)  An  understanding  of  the  factors  controlling  the  behavior  of 
detached  flows  is  a prerequisite  to  being  able  to  predict  it.  The 
studies  of  Cp  and  flow  visualization  of  diffusers  by  the  Stanford  HTTM 
grouo  over  the  last  15  years  have  greatly  increased  the  understanding  of 
the  qualitative  features  of  these  flows.  These  studies  nov  need  to  be 
extended  to  include  detailed  quantitative  flowfield  information,  such  as 
the  turbulence  field,  intermittency  and  skin-friction  along  the  walls. 
Because  of  the  complicated  nature  of  the  detached  flow  regions,  these 
measurements  will  not  be  easy,  and  new  measurement  techniques  such  as 
laser  Doppler  anemometers,  etc.,  may  have  to  be  developed. 

(b)  The  currently  used  eddy  viscosity  concept  is  a gross  approxima- 
tion to  the  actual  flow.  When  new  data  become  available,  scaling  laws 
relating  the  shear  stresses  to  the  turbulent  kinetic  energy  or  entrain- 
ment will  permit  improved  calculation  methods  to  be  developed. 

(c)  The  current  method  can  be  extended  quite  readily  to  the  1-D 
core  axisymmetric  case  for  both  incompressible  and  compressible  diffu- 
sers. The  next  step  is  the  case  with  the  incompressible  inviscid  core 
calculated  from  a solution  of  Laplace's  equation  in  the  axisymmetric 
effective  flow  channel.  The  corresponding  compressible  case  must  await 
the  development  of  a fast  algorithm  for  compressible  potential  flow. 

(d)  An  alternative  approach  to  the  iterative  matching  procedure 
between  the  boundary  layer  and  the  core,  and  the  inherent  convergence 
problems  thereof,  is  to  couple  both  regions  into  one  large  domain  and 
solve  the  whole  flowfield  as  an  elliptic  problem.  The  equations  for 
such  a scheme  have  been  developed,  but  no  solution  has  been  attempted. 

The  approach  looks  promising. 
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Straight-walled  diffuser  flow-regime  chart  of  Fox  and 
Kline  [1]. 
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Fig.  2.  Ability  of  Eqn.  (2-15)  to  represent  attached  boundary  layer 
velocity  profile?.  Data  are  from  Strickland-Simpson  [32], 
station  88.2  inch. 
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Fig.  3a.  Ability  of  Eqn.  (2-15)  to  represent  detached  boundary  layer 
velocity  profiles.  Data  are  from  Strickland-Simpson  [32], 
station  157.1  inch. 
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Fig.  3b.  Ability  of  Eqn.  (2-15)  to  represent  detached  boundary  layer 
velocity  profiles.  Data  are  from  Tani's  [38]  flow  over  a 
backward-facing  step. 
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Fig.  7.  Entrainment  Equation  Summary 
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Fig.  8.  Comparison  of  Tmax/pU°°  and  entrainment  rate  data  with  that 
obtained  from  Eqn.  (2-40).  Data  are  from  Str ickland-Simpson 
[32].  Intermittent  detachment  is  at  x = 127  in  and  x = 0 
at  x = 132  in. 
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Fig.  9a.  Effect  of  lag  parameter  Aa  on  Bradshaw-Ferr iss  (2400)  relaxing 
flow  (a  * -.255  -*•  0).  Prescribed  pressure  gradient  calculation 
Mean  boun  lary  layer  parameters. 
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Fig.  9b.  Effect  of  lag  parameter  Xa  on  Bradshaw-Ferriss  (2400)  relaxing 
flow  (a  = -.255  -*•  0).  Presecribed  pressure  gradient  calculation. 
Skin  friction  and  entrainment. 


Fig.  12.  Results  — Bradshaw-Ferriss  (2600)  equilibrium  flow  (a 
Prescribed  pressure  gradient  calculation. 
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Fig.  13.  Results  — Clauser's  equilibrium  flow  (2200)  in  mild  positive 
pressure  gradient.  Prescribed  pressure  gradient  calculation. 
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Fig.  16.  Results  — Perry  diffuser  flow  (2900)  showing  comparison 
between  prescribed  pressure  gradient  and  the  1-D  core 
diffuser  calculation. 
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Strickland-Slrapson  flow  (lower  wall)  as  calculated  with 
prescribed  pressure  gradient.  (a)  x = °°,  (b)  x = 

I indicates  location  of  intermittent  detachment. 
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Fig.  19.  Results  — So-Mellor's  [46]  convex  wall  boundary  layer  as 
calculated  with  prescribed  pressure  gradient. 
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Effect  of  lag  parameter  Aa  on  the  predictions  of  the  unstalled 
diffuser  flow  of  Carlson  et  al.  [27].  Calculation  is  1-D  core 
with  symmetric  boundary  layers. 


Effect  of  lag  parameter  Xj  on  the  calculations  of  a diffuser 
in  transitory  stall  as  measured  by  Carlson  et  al.  [27]. 
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Predicted  variation  of  boundary  layer  quantities  H,  6 and 
Cf/2  along  the  walls  of  a diffuser  operating  in  the  transitory 
stall  regime.  Diffuser  is  same  as  from  Carlson  et  al.  [27], 

Run  62430,  N/V^  « 6,  AR  = 2.4,  B1  - .030.  No  data  are 
available  for  comparison. 
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Predicted  exit  conditions  for  N/W^  ”6,  Bi  - .030  diffuser 
family.  Only  one  data  point  is  available  for  comparison. 

Data  are  from  Carlson  et  al.  [27],  Run  62430. 
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Fig.  30.  Strickland-Simpson  [32]  flow,  comparing  data  with  predictions.  Full  2-D  core 
solution  in  inviscid  core  from  x = 8.11  ft  to  exit  of  diffuser.  Prescribed 
pressure  gradient  calculation  from  inlet  to  x = 8.11  ft,  at  which  point 
(marked  A)  the  upper  boundary  layer  was  removed.  No  lag  equation  in  reversed 
flow  region. 
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Fig.  31.  Strickland-Simpson  [32]  flow.  Same  run  as  Fig.  30. 
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USER'S  GUIDE  TO  PROGRAM  TSTALL 


UG1 . INTRODUCTION 

TSTALL  is  a FORTRAN  program  that  performs  four  types  of  computa- 
tions. 

(a)  Turbulent  boundary  layer  development  with  prescribed  pressure 

(b)  Calculation  of  diffusers  operating  in  the  unstalled  and  transitory 

stall  regimes  with  20/20  = 1.2,  assuming  one-dimensional  core 

and  symmetric  b.l.'s. 

(c)  Calculation  of  the  same  diffusers  as  in  (b)  but  with  a two- 
dimensional  core  assumption. 

(d)  Solution  of  Laplace's  equation  using  a boundary  integral  method. 

Access  to  the  subroutines  for  each  calculation  is  done  in  the  MAIN 
routine  according  to  inputted  keywords  for  geometry,  problem,  and  core 
types.  The  options  available  are: 


Geometry 

(GEOMT) 


type 


Standard  symmetric  diffuser  (STDD) , where  the  co- 
ordinates defining  the  geometry  are  internally 
generated. 

Non-standard  diffuser  (NSTD) , where  the  user  enters 
the  geometry. 

Standard  asymmetric  diffuser  (HALF),  where  again 
the  geometry  is  internally  generated. 


Problem  type 
(PROBT) 


Calculate  the  TBL,  given  the  pressure  gradient 
(TBLP) . 

Two-dimensional  inviscid  calculation  only  (NOBL). 
One-  or  two-dimensional  diffuser  calculation  (DIFF). 


One-dimensional  core  (ONED). 

Core  type  

CORET)  ~ Two-dimensional  core  (TWOD) . 


Flowchart  for  the  MAIN  routine  is  shown  in  the  figure  below.  Names 
of  called  subroutines  are  marked  above  the  relevant  boxes. 


U1 


t 


UG2.  GENERAL  CONSIDERATIONS 

The  input  data  and  output  for  each  type  of  computation  is  presented 
in  the  next  few  sections. 

All  input  format  field  widths  are  E10.0  for  real  variables  and  110 
for  integers.  Integer  inputs  must  of  course  be  right-justified. 

Multiple  runs  can  be  made  by  stacking  the  corresponding  input  cards. 
The  program  recognizes  two  blank  cards  as  the  end  of  input  data. 

Either  FPS  or  MKS  units  may  be  employed.  Data  entered  in  either 
system  will  result  in  output  of  the  same  type.  For  example,  SWU(m) , 
VISC0S(m2/sec) , etc.,  will  give  output  X(m) , DSTAR(m) , UB(m/sec),  etc. 

The  diffuser  is  divided  into  N segments  which  are  numbered  in- 
creasing in  the  counterclockwise  direction,  as  shown  in  Fig.  U6.  The 
node  number  0 is  assigned  to  the  lower  wall  inlet  location.  Since  zero 
subscripts  are  not  allowed  in  FORTRAN,  the  zeroth  node  is  reassigned  the 
value  of  N,  and  is  equal  to  36  for  the  sample  case  shown  in  Fig.  U6. 

The  program  does  this  automatically  when  GEOMT  is  set  equal  to  STDD  or 
HALF.  When  using  the  nonstandard  option  (NSTD),  the  user  must  adhere 
to  this  numbering  system. 

Internally  generated  segments  allow  for  a non-constant  node  spacing 
to  allow  for  greater  resolution  in  this  region.  The  points  are  spaced 
according  to  a geometric  progression,  the  spaces  increasing  in  both  direc- 
tions away  from  the  throat. 

Default  values  have  been  used  wherever  possible,  and  may  be  used  as 
indicated  on  the  input  card  image  data  outlined  next.  A blank  or  zero 
entry  will  result  in  the  use  of  the  default  value. 

Diffusers  that  are  to  be  calculated  using  this  program  must  meet 
the  following  requirements: 

(a)  Aspect  ratio  AS  •>  5. 

(b)  Straight  walls,  since  the  b.l.  cannot  handle  curvature. 

(c)  Inlet  blockage,  Bj  £ .050,  and  flow  at  inlet  must  be  turbulent 
and  1-D. 
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UG3.  B.L.  CALCULATION  WITH  PRESCRIBED  PRESSURE  GRADIENT 


The  input  data  may  be  conveniently  entered  using  the  template 
shown  in  Fig.  Ul.  A description  of  the  input  parameters  on  each  card 
follows,  with  the  format  information  in  parentheses  at  the  end. 


Card  1.  Title  for  user  identification  (A(80)). 


Card  2. 


Card  3. 


Card  4. 


Card  5. 


Card  6. 


Card  7 . 


Card  8. 


Keywords  specifying  geometry  (NSTD)  and  problem  (TBLP),  as 
shown  (3(6X,A4)). 

Number  of  stations  along  the  flow  for  which  velocity  data  will 
be  inputted  (110). 

Starting  values  of  XX,DELST,H  and  the  kinematic  viscosity 
VISCOS  (4E10.5). 

Interval  between  b.l.  printouts,  IPR  (recommended*l) , location 
of  3-D  source  XC  (default=lE5  if  left  blank). 

Repeated  values  of  XX  station  location  SWU  and  the  correspond- 
ing velocity  VIU,  until  all  data  are  exhausted  ( 8E10.Q ). 


Sample  input  for  Bradshaw-Ferr iss  relaxing  flow  (2400)  is  shown  in 
Fig.  U 2.  The  corresponding  output  is  presented  in  Fig.  U3  and  plotted 
in  Figs.  9a  and  9b. 
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CALCULATION  WITH  PRESCRIBED  PRESSURE  GRADIENT 


Template  for  turbulent  boundary  layer  calculations  with  prescribed 
pressure  gradient. 
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UG4 . STANDARD  DIFFUSER  CALCULATION  WITH  1-D  OR  2-D  CORE 


The  template  shown  in  Fig.  U4  details  input  data  for  diffuser  cal- 
culations with  either  a one-  or  two-dimensional  inviscid  core.  Refer  to 

Figs.  U5  and  U6  for  standard  diffuser  geometry  and  nomenclature. 

Card  1.  Title  for  user  identification  (A(80)). 

Card  2.  Keywords  specifying  geometry  (HALF  or  STDD) , problem  (NOBL  or 
DIFF),  and  core  type  (ONED  or  TWOD) . (3(6X,A4)). 

Card  3.  XI,  RC1 , N,  RC2 , X2,  W1  (Fig.  U5) , TWOTHD  (26  in  degrees  — for 
asymmetric  unit,  enter  twice  the  angle  of  the  diverging  wall), 
aspect  ratio  AS  (default  AS-8) . (8E10.5). 

Card  4.  Number  of  segments  in  inlet  (Nl) , throat  curve  (NCI),  diffusing 
section  (N2) , exit  curve  (NC2) , tailpipe  (N3).  ND1  and  ND2  are 
fractions  of  the  inlet  and  diffusing  sections  where  node  nearest 
the  throat  is  to  be  located  (default  ND1=5*N1,  ND2=5*N2) . (7110). 

k 

Card  5.  Inlet  blockage  Bl=26  /Wl,  inlet  velocity  Uco(UI),  kinematic  vis- 
cosity VISCOS,  location  of  source  or  sink  for  3-D  correction  XC 
(default  XC=1E5).  (4E10.5). 

* 

Card  6.  Inlet  b.l.  parameters,  lower  wall  H and  6 (HS.DELSTS),  and 

upper  wall  (HU,  DELSTU) . Leaving  HU  and  DELSTU  blank  implicitly 
sets  them  equal  to  the  lower  wall  values.  (4E10.0). 

Card  7.  Boundary  layer  print  interval  IPR  (recommended=2) , type  of  print- 
out for  interior  points  in  the  inviscid  core  NORMPR  (for  NOBL 
option  only,  =-l,+l,0  for  normalized,  dimensional  or  both), 

1TMAX  the  maximum  allowable  iterations  for  2-D  core  diffuser 
calculation  (recommended  8 to  10),  CPEROR  is  the  convergence 
criterion  (recommended  .025).  ( 3110 , E10 . 0) . 


Sample  input  for  the  1-D  core  diffuser  of  Carlson  and  Johnston  [27] 
AR=2.4,  N/W1=6,  Bl=.025,  is  shown  in  Fig.  U7.  The  corresponding  output 
is  shown  in  Fig.  U8  and  plotted  in  Fig.  21. 

For  inviscid  calculations  (PR0BT= 'NOBL' ) , additional  input  data  are 
needed  for  determination  of  the  values  of  the  complex  function  and  its 
derivatives  at  interior  points. 

Card  8.  LINES  is  the  number  of  lines  along  which  interior  point  values 
are  to  be  computed  (default  0).  (110). 
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Card  9.  Enter  one  card  for  each  line  along  which  interior  point  valu 
are  desired.  XI, Y1  (coordinates  of  start  of  line),  X2,Y2  (en. 
of  line),  NSEGS  (number  of  segments  that  each  line  is  divided 
into).  Do  not  place  any  interior  point  on  a node  location. 
(4E10. 0,110) . 


Card  10. 


Card  11. 


Sample  input  for  an  inviscid  solution  of  the  Carlson  et  al.  diffuser 
is  shown  in  Fig.  U9  and  the  corresponding  output  in  Fig.  U10. 


STANDARD  DIFFUSER  CALCULATION  WITH  1-D  OR  2-D  POTENTIAL  CORE 


Fig.  U4.  Template  for  standard  diffuser  calculations.  Cards  8 through  the  end 
are  for  interior  point  computations  only  and  are  not  required  for  dif- 


Segment  locations  for  a typical  36-node  standard  diffuser. 


**CARLSON  RUN  62430  ,N/W1=6.0,  Bl=0.025,  Wl=.2500,  H=1.41,  2THETA= 13 . 3 09 
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UG5.  NONSTANDARD  DIFFUSER  CALCULATION  WITH  1-D  OR  2-D  CORE 

For  nonstandard  diffusers,  node  points  describing  the  geometry  have 
to  be  entered  by  the  user.  The  location  of  the  axes  and  the  numbering 
system  must  be  the  same  as  that  for  the  standard  case.  Fig.  U6.  The  tem- 
plate, Fig.  Ull,  describes  the  input  data. 

Card  1.  Title  for  user  identification  (A(80)). 

Card  2.  Keyword  specifying  geometry  (NSTD),  problem  (NOBL  or  DIFF) . and 
core  type  (ONED  or  TWOD) . (3(6X,A4)). 

Card  3.  Number  of  segments,  total  N,  lower  wall  NR,  upper  wall  NL.  For 
a 2-D  core  calculation,  NR  should  be  = NL.  N=NR+NL+2.  (3110). 

Card  4.  Enter  (XW,YW)  coordinates  of  segment  end  points,  and  ALW  the 

angle  between  the  segment  and  the  positive  X direction  (CCW  posi- 
tive CW  negative).  Enter  one  card  for  each  node,  beginning  with 
node  1 and  ending  with  node  N (which  is  really  node  zero).  (3E10.0) 

Card  5 

Card  6 


Card  (4+N) . 

Card  (5+N).  Inlet  width  Wl,  divergence  angle  TWOTHD  (degrees),  aspect 
ratio  AS  (default  AS=8) . (3E10.0). 

Card  (6+N).  Inlet  blockage  Bl,  inlet  velocity  UI , kinematic  viscosity 
VISCOS,  location  of  3-D  source  or  sink  XC  (default  1E5). 
(4E10.0). 

Card  (7+N).  B.l.  print  interval  IPR  (recommended=2) , type  of  printout  for 
interior  points  in  the  Inviscid  core  NORMPR  (for  NOBL  option 
only,  =-l,+l,0  for  normalized,  dimensional,  or  both),  ITMAX 
the  maximum  number  of  iterations  allowed  for  2-D  core  calcula- 
tion (recommended  8 to  10).,  CPEROR  is  the  convergence  criter- 
ion (recommended=.025) . (3I10.E10.0) . 

* 

Card  (8+N).  Inlet  b.l.  parameters,  lower  wall  H and  6 (HS,DELSTS),  and 
upper  wall  (HU,  DELSTU) . Leaving  HU  and  DELSTII  blank  im- 
plicitly sets  them  equal  to  the  lower  wall  values.  (4E10.0). 

Fig.  U12  is  sample  input  data  for  a 2-D  core  calculation  of  Strick- 
land-Simpson ' s flow.  Fig.  U13  is  the  first  and  last  few  pages  of  output, 
which  are  plotted  in  Figs.  30  and  31. 

If  only  an  inviscid  solution  is  desired,  cards  8 through  the  end  of 
the  standard  diffuser  input  should  be  added  to  the  end  of  this  deck. 
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NONSTANDARD  DIFFUSER  CALCULATION  WITH  1-D  OR  2-D  POTENTIAL  CORE 


Template  for  nonstandard  diffuser  calculations 
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//STAND  JOB  'J1S«D1,511' .‘S-GHOSE' ,CLASS=E,TIHK=(,90) 
//  FTEC  E0RTC1,PAPH.F0RT=' 0PT=  2* 

//EOPT.SYSIN  DD  * 


C MAIN  PO!I  TI  HE  TO  CALC'JI.ATF  THE  P FR  FOR  HA  HC  E OF  1-D  AN  I 2-D 

C DIPEOSBRS,  OPERATING  IN  THF  UNSTALLFD  AND  TRANSITORY  STALL 

C REGIMES.  THE  SCHEME  OSES  A NEK  TO  RBIILENT  BOUNDARY  LAYER 

r °PED TCTION  METHOD,  TOGETHER  WITH  SIMULTANEOUS  ITERATION 

C RETHEFN  THE  CORE  AND  THE  BL.  AN  ADDITIONAL  INVISCID  BATCHING 

C TECH  NT  OWE  IS  USED  TO  HATCH  THE  2-D  FLOWFIELD  KITH  TEE 

C THE  STPEAHTUBE  ENVELOPING  THE  B L. 

C WRITTEN  PY  SANJOY  GHOSE,  MECHANICAL  ENGG  DEPT,  STANFORD  ONIV. 

C LAST  REVISION  HADE  ON  NOV  1,1976. 

C 


INTEGEP  HEAD  (20)  , GFOH (4)  , CORE (4) , PROB  (4 ) , GEOH  T, COR  ET  ,PROBT# I C 1, 

$ I D2,  IDS 

COHHON/OD E 1 S/JSTPTS ,JENDS, NDTM,SW( 90)  , HI (90)  ,DWI (90)  ,DDWI(90) , 

*DS  (90) 

COMMON /N  ST  D/TE1, 102, ID  1, NST, S HT  (90 ) 

DATA  GEOH/ * STDD ' , ' • , • NSTD • , • HA LF' /, COP E/' ON  ED* , * THO D' , 

♦'  TXXX  '/.PROD/'TBI.P'  , • NOBL  ' , • DIFF*  , ' •/ 

NDTH=9D 

IDO  READ (S, 900, END=B0O)  (H E AD  ( J)  , J= 1 , 20 ) 

WRITE  (6, 9 10)  (HEAD(J)  ,J  = 1,2D) 

P^AD  (6,97  0) GEOHT ,P ROOT, COPET 

W R IT  E ( 6,  910)  G EOMT  , PRO  BT,  COR  ET 

rpl  = 0 

T D2  = 0 

TD1=0 

NST  = 0 

no  1 OS  J = 1 ,4 

I F (GEOH (J)  . E0. G^OHT)  ID1=J 
TF  (PROO(J)  . EO  • PROPS’)  TD?=  J 
IOC  IP  (CORE  (J)  .F0.COPFT)  I 03= J 

TP(MAX0  (ID1,ID2,ID1)  ,GT.0)GO  to  1 1 0 

HRTTE  (6,940) 

STOn 

110  IF  (ID1  .so.l)  GO  TO  112 
I’’  CALL  STAND 
GO  to  1?n 

112  I0,  (TD2  .EO  . 1)  GO  TO  120 
N ST=  1 

CALL  N ST  AN  D 

’in  rn  TO  (l?i, 122, 179, 12’) ,ID2 
121  CAM.  PSTEST 
GO  TO  10O 
12"  CALI.  INVCTD 
GO  Tp  109 

1 2 1 TFMDl.EO.I)  C ALT  DIFE1D 
T" (TD1.«F. 1) C ALL  DI Fw? D 
GO  TO  1 o o 
900  w PTTp  (6,  9 SO) 

STOP 

900  FORMAT  (20 AU)  .1 

91"  ®OR«AT  ('  1 • ,20  A4//) 

9 ?0  c O PM  A t f 7 ( 6X  , A 4)  ) 

910  FOP" AT  ( • GFOHTTPYs  ','9,'  pROPLEH  type*  *,A4, 

**  cort  velocity  profile*  »,au//) 

040  "OPMATf*  ♦••"CANNOT  RFCOGHTEF  PPOPLFH  TY  PF- -CHECK  CAPE  • 2 ****//) 
qeo  CO- 1 »-(  M *♦*♦♦• r N O Of  PROG  PA M *•••••• ) 

ph  n 


U14 


s tii so ttt t n 17  adans  ( ny  ,neo,  dnaie,  irungf) 


r USES  UmH  OP  P FR  A PAH  S-NOULTON  PF FDTCTOP-CORPFCTOR  H FT  HOD 

r SO  S OI.  V E A S ET  OP  FIRST  ORDER  CDF'S  EXPRFSSFD  IN  THE  PORN 

C Y'=P  (X.Y....  .)  . OSES  A 4TH  ORDER  PUNGE-KUTTA  NFTfcOD  FOR 

C STARTING  {AND  RESTARTING) . CALL  TO  THIS  ROUTINE  REURNS  VALUES 

r op  THF  FUNCTION  A^  XfDX,  GIVEN  VALUES  AT  X. 

r RATE  NATRIX  CONTAINS  DERIVATIVES  POR  THE  LAST  4 STFPS. 

C TH p ROW  P ATE  (1, J) , J= 1 ,NEO  HAS  VALUES  FOP  STEP  N, 

C POW  R ATE  { 2, J)  FOP  STEP  (N-I).PTC.  VA  LS (J) , J = 1 , NFQ  ABE 

C VALUES  OP  THE  VARIABLES. 

C SET  IPUNGE=1  IN  THE  CALLING  POUTINF  TO  PROVIDE  STARTING 

C VALUES  VIA  CALL  TO  PKS4.  IRUNGE=S  CAUSES  ADAHS-HOULTON 

C «th  ORDER  PREDICTOR-CORRECTOR  M FT  HOD  TO  BE  INVOKED. 

C 


EXTPRNAL  PNANE 

PEAL  VALP  (8)  , VALC  (8)  .RATEP  (8) 

CONNON/ADAH1/X  , VALS  (4)  , RATES  (4)  .RATE (4, 8) 

rONNON/ODEIS/JSTRTS,  JENDS,  NDIH.SH  (90)  , HI  (90)  ,DWI(9P)  ,DDHI(90)  , 
*ns (90) 

DATA  STFPER/1  . E-V 
IF  (NE0.LE.8)  GO  TO  10(1 
WRITE  (6,900)  NEC 
STOP 

IOC  CONTINUE 

RRRT  OW=STEPER/S.O 
IK  IP  (TRHNGE.EO.S)GO  TO  200 
CALL  RKS4  (DX.NEO.DNANE) 

DO  140  J = 1 , N FQ 
140  RATF (TPUNGE.J) =RATES (J) 

TPUNGP=TRUNGE*1 

RETURN 

C 

C START  ADAHS-BASHFORTH  PREDICTOR  ROUTINE 

200  X =X ♦ DX 

DH=DX/24 . 0 
DO  210  J=1,NEQ 

V AIP  (J)  = V ALS  (J)  ♦0H*(55.0*FATE  (1.J)  -59  ,0*R  ATE  (2  , .1)  ♦ 37 . 0*RATE  ( 3,  J) 
S-9.0FRATE  (4  , J)  ) 

230  CONTINUE 

L i c 

CALL  DNANB(X.VALP.RATEP) 

C 

C BEGIN  ADAHS-NOULTON  CORRECTOR  ITERATION. 


in 


RATES(J)«RATBP(J) 

105  RATE (1,J) =RATBP(J) 

I R0RGB*5 

IF(DBRR.GT.ERBLOW) RETORR 
DX*2. 0*DX 
IRURGE=1 
RETORR 

OHABLE  TO  CONVERGE  IR  2 ITERATTORS  OF  THE  CORRECTOR.  COT 

STEPSIZE  TR  HALF  AND  RESTART. 

310  X -X-DX 
DX=0. 5*DX 
I RUNGE=1 
GO  TO  110 

900  FORHAT  ( * NOHBER  OF  EQOATIOHS  EXCEEDS  ABRAT  BOUNDS,  AS  NEQ=  ',15, 
S'  INCREASE  DI-HEN  SION S OF  SCRATCH  ARRAYS  IR  ADAHS  AND  RKS 4 •) 

END 

SUBROUTINE  DIFF1D 

C ROUTINE  TO  TEST  1-D  CORE  DIFFUSERS  IN  SIHOLTANEOUS  ITERATION. 

COHHON/DER2/DELT, UB,UT,DI,VT,VB,UDUI,TAUH,H,THETA,DELST,CPD2, 
SVISCOS.NBL 

COHHON/ODE1S/.TSTRTS,JENnS,NDIH,SW(90) ,W I (90)  , DW I (90) ,DD»I  (90) , 
$ns  (90) 

CONNOH/SPLYN/XX,VT,DWT,DDWT,ISFTDP,KHID 
COHHON/TE  HP1/XCHX , IN  ALLV 
XX=0. 0 
NBL=2 

CALL  TBT.SI(O) 

WRITE  (fi,93r>) 

930  FORH AT (' 0 **********  END  OF  ROUTINE  DIFF1D********') 

RETURN 

*•90 

SUBROUTINE  F ACTOR ( A , W , I PIVOT , D, N , IFLA G) 

C 

C THIS  SUBROUTINE  PERFORHS  A L-U  DECCH POSITION  ON  TEE 

C GIVEN  HATRIX  A (N , N) , AND  RETURNS  THE  HATPIX  V.  IPIVCT 

C . IS  A VECTOR  CONTAINING  THE  PIVOTING  ORDER. 

C 

SEAL  A (N,  N)  ,W  (If,  N)  ,D(N) 

INTEGER  IPTVOT(N) 

TFLAG=1 

C INITIALIZE  B , IP! VOT , D 

DO  10  1=1 ,N 
I PIVOT  (I)  = I 
ROWH AX=0 . 0 
DO  9 J = 1,N 
H f I,  J)  = A ( I , J) 

ROBHA  X = AH  AX1  (ROWH AY , A BS ( W ( I, J) ) ) 

9 CONTINUE 

IF (ROBHAX. EQ. 0.0) GO  TO  999 
D (I)  = ROWH  A X 
10  CONTINUE 

C GAUSS  ^LIHINATION  WITH  SCALED  PARTIAL  PIVOTING 

NH 1 = N - 1 

IF  (NH 1 .E0.0) PETURN 
no  20  |r=1  , NHl 
3 = K 

K PI =K ♦ 1 
Ip=TprvoT  (F) 

COLH  A X =A  P S (W  (TP,  K)  ) /D(IP) 
no  11  t=*P1,N 
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TP=TPI  VAT  (T) 

AWTKOV=AFS(W  ( TP,  K ) ) /D  ( TP  ) 

TPfAWTfnV.LF.CnLMA'f)  GO  TO  11 
COT.NAX  = AWIKOV 
J = T 

11  CONTINUE 

TP(COLMAX. Fo. 0.0)  RO  TO  999 
r 
c 

T PK=I PITOT (J) 

T PIVOT  (J) =l?TVOT (K) 

TOTVOT  (K)  = IPK 
DP  20  T=KP1,N 
T P=T  P I VOT  (T) 

W (TP,  K) =W  (TP,K)/W (IPK,K) 

9 ATT 0 = -W  (TP,K) 

90  20  <1=FP1,N 

W (TP,.7)=PATT0*W(IPK,J)  + W(TP,J) 

20  O0VTTNUp 

TP  (W  (TP,  9)  .B0.O.0)  00  TO  999 
p pTrip  n 
C 

C 9FT  IFT.  AG  = 2 TO  INDICATE  INAPTLITY  TO  FACTORIZE  MATRIX. 

C 

999  TFLAG=2 

W PT?  F (6,99991 

9909  PORIATf  1 HO, • ****ONABLE  TO  COMPLETE  L-0  DECOMPOSITION  OF  HATRIX*) 

STOP 

END 

SUBROOTTNE  PpSL 
9FAL*9  A (96,90) 

COMPLEX  C (91) 

COMPLEX  CMPLX.ICMPLX 

RPAL  LNV(90)  ,VEL(9  0)  ,70(120)  , TO  (12  0) 

COHMON/FPCV  AL/7C (90 ) ,YC(90) ,AL(90)  ,LNV 

C OMMON/G  EOM 1/X  D,  XL , T H, WIDTH, X 1 ,B 1, SINTH ,COSTH , T BOTH R ,TWOT HI , 
SSTNTH1 ,rOSTH1 ,AS,WH,XSP,X2,XHAX,XDE 
COMHO  V /G  EOH2/N , NR , NL,N0,NM1, NLC , NPC 
COMA  "«/TN  IT  AL  /DEL  ST  1, HI,  U II, I PR  1 
COHMO  S/PF SL1 /C, A , XO , TO 


r N = NO  h r FT?  OF  SEGMENTS  (PROH  STAND  OR  NSTAND)  . XC(J),TC(J),  AND 

C A I.  (.1)  , ARE  VALDES  AT  THE  EDGE  OF  THE  EFC  AS  PASSED  VIA 

C COMM  ON/E  FCV  AL/. 

C NA  ROW  = N0  OF  ROWS  OF  A,  NACOL*NO  OF  COLS  OF  A. 

C 


"Al'A  NAROW,NACOL,  VI  NORM/86, 87, 1 .0/ 
ND=N-2 
N0P1=N0*1 
V SCA  LE=DT 1 


C COMPOTE  THF  BOUNDARY  COORDTNATFS  IN  THE  COMPLEX  P LA  KE. 

DO  200  .7=1  ,N 

200  C (J)  =CMPLX  (XC  (J)  ,YC(J)  ) /XL 
C ASSIGN  STARTING  VALDES 


LNV  (N)  =0.0 
LNV  (NN1 ) =0.0 
CALL  EEPSEG 

CALL  GJR  (A, NO,  1.  E- 1 0,  N AROW,  N ACOL) 

A (NO  * M0P1 ) IS  MATRIX  OF  COEFFS  WITH  B VECTOR  STORED  IN  LAST  COL, 

ANSWER  LNV  VECTOR  RETnRNED  IN  LAST  COL,  NUP1. 

DO  72 S .7  = 1 , NO 
72S  LNV(.I)  = A ( J,  NOP  1) 
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DO  810  J=1,N 
810  VBL  (J)*EXP(LWV(J)  ) 

VEWORH=  (TEL (NLC)  ♦TEL(BRC)  ) /2. 0 
WRITE (6,940)XL,VSCALE,VINORH,VEWORH 
WRITE  (6, 9 4 2) 

DO  883  J=1,N 
CP*1. 0-VEL  (J)  **2 

883  WRITE  (6,94  3)  J,C(J) ,AL(J)  ,LNV(J) ,VBL(J),CP 
RETORB 

940  F0RHAT(1H1,46X,' LENGTH  SCALE ', 1 4X, •=', 1 PE14 .5/1 H ,46 X ,' VELOCITY  SC 
SALE'  ,12T,  •=•  ,1PE14.5,/1H,46X,'NORHALIZED  INLET  VELOCITY *1PE 
S14.5/1H,46I,'NORHALIZED  EXIT  TELOCITI', • = • , IP  E 14 . 5/  IfiO,  5 3X, 
5'NORHALIZED  SOLUTION') 

942  FORHAT(1H0,'i',9X,'XC',12X,' YC ,12X,'  ALPH A' ,7  X , 'LN ( V It)  ',6X, 

F,  • TELOCIT  T',5X,'CP') 

943  F0RHAT(1H  ,I3,6E14.6) 

END 

SOBROOTINE  RKS4 (H,» EQ , DR  AH E) 

SOLTE  A SET  OF  FIRST  ORDER  ODE'S,  T • = F (X , T ( 1) , Y ( 2) , . . I (NEO) ) 

USING  FOURTH  ORDER  RUNGE-  KUTTA  SCHEHE  WITH  FIXED  STBPSIZE  H. 
BETURNS  TALOES  OF  VECTOR  I AT  X*DX  GIVEN  VALOES  AT  X. 


REAL  KV(8,4)  ,SPAN(4)  ,10(4) 

COHHON/AD AH  1/X, T (4)  ,F(4)  ,RATE(4,8) 

EXTERNAL  DNAHE 

DATA  SPAN/O.  5,0.  5,  1.0,  1. 0/ 

XO=X 

DO  200  J=  1, N EQ 
?0o  TO  (J)  =T(J) 

DO  400  1=  1,4 
CALL  DNAHE  (X,  I,  F) 

DO  300  J = 1 , N EQ 
300  KV  (I,  J)  =H*F  (J) 

DO  350  J= 1,NEQ 

350  T (J)  =10  (J)  ♦SPAN  (I)  *KV  ( I,  J) 

400  X=XO*H*SPAN (I) 

DO  500  1=1, NEQ 

50  0 T (I)  =T0  (T)  ♦ (K  V ( 1 ,1)  + KV  (4,1)  ♦2.0*(KV  (2,1)  ♦ KV(3,I) ) )/6.0 
RETURN 
END 

SOBROOTINE  SUBST  ( H ,B ,X  , IPIVOT , N) 

PERFORH  BACK  AND  FORWARD  SUBSTITUTION  TO  CALCULATE 

THE  UNKNOWN  VECTOR  X,  AS  THE  SOLUTION  OF  A*X=B. 

REAL  * (N,  N)  ,B  (N)  ,X(N)  ,SUH 
INTEGER  IPIVOT  (N) 

I F (N. GT. 1 )GO  TO  30 
X (1)  =B  (1)  /W(1  ,1) 

RETURN 

30  I P=T PT  VOT  (1 ) 

X (1)  =B(TP) 

DO  50  K=2 , N 
I P=I PIVOT  (K) 

KH 1 =K- 1 
S UH=0 . 0 
DO  40  J=1 , KH 1 
40  SUH=W  (TP,  .1)  *X  (J)  ^SUN 
50  X (K1  =8  (TP)  -SUH 


X(N)=X(N)/W(IP,N) 
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on  60  J=K  PI , N 
60  SUM=H  (TP,  .1)  *X  (J)  +SUM 
IB  X (K)  = (X  ( K ) -SUN)  /H  (T  P,  K) 

RETURN 

END 

SUBROUTINE  T RID  AG ( I F,L , X , B ,C ,D , V, 0 , HDT H) 

r SUBROUTINE  FOR  SOLVING  A SYSTEM  OF  LINEAR  SIMULTANEOUS 

r FOU  AT  ION  S HAVING  A TRIDIAGONAL  COEFFICIENT  MATRIX. 

C THE  EQUATIONS  ARE  NUMBERED  FROM  IF  THROUGH  L,  AND  THEIR 

C SUB-DTAGONAI,,  DIAGONAL,  AND  S UPE R-DI AGON AL  COEFFICIENTS 

C ARE  STOPED  IN  THE  ARRAYS  A, B,  AND  C.  THE  COMPUTED 

C SOLUTION  VECTOR  V(IF)...V(t)  IS  STORED  IN  THE  ARRAY  V. 

PFAL  A (M)  , B (M)  ,C  (M)  ,D(M)  ,V(NDTN)  ,BETA  (10  1)  , CAMN  A ( 10  1) 

C ...COMPUTE  INTERMEDIATE  ARRAYS  BETA  AND  GAMMA... 

BFTA  (IF)  = B (IF) 

GAMMA  (IF)  = D (IF)  /BETA  (IF) 

TFP1  = I F*  1 
BO  1 I=IF  P 1 , L 

BETA  (I)  = B (I) -A  (T)  *C(I-1) /BETA  (1-1) 

1 GAMMA(I)  = (D  (II  - A (I)  *GAMMA  (1-1)  )/BETA  (I) 

C ...COMPUTE  PINAL  SOLUTION  VECTOR  V... 

V (L)  = GAMMA  (L) 

LAST  = L-TF 
DO  2 K=1 , LAST 
I = L-K 

2 V(I)  = GAM  HA  (I)-C(T)  *V(I*1)/BFTA(I) 

RETURN 

END 

SUBROUTINE  CHANGE 

COMMON/GEOM2/N,NR,NL,NH,NM1, NLC, NRC 
CONNON/TENP1/XCMX,IWALLV 
COMMON/CON/SV AL (RB ) ,YVAL(90) 

COHMON/ODE1S/JSTRTS,JENDS,NDIM,SH(90)  ,WI (90) , DHI  (90)  ,DDHI  (90)  , 
*DS  (99) 

COMHO N/O  DEI U/JSTR  TU,JENDU,SNU (90) ,HIU (90) 

COMMON/ODE2U/JTBL  U,STBLU  (90)  ,DSTARU(90)  ,UI1D0(90)  ,DELTU(90)  , 
*T!T2PU  (9B) 

COHMON/O DE2S/JTBLS , STBLS  (90)  , DST  ARS  (90)  ,0I1DS  (90)  ,DEITS(90)  , 
SUT2DS  (9B) 

COHMON/SPLYN/XINT, PINT , FP TNT , FPP INT,I SETU P, KHI D 
RFAL  HELP  (90) 

NMNLC=N~NLC 

C SET  UP  THE  SPLINE  COEFFICIENTS  FOR  SVAL,  YVAL. 

XTNT=o.n 

ISETnp=0 
M HID=  2 

r INTPRPOLATF  FOR  THE  VALUES  OF  YVAL  AT  THE  HALL  LOCATIONS  SVAL. 

IVIHALI.V.  EQ.  1)  GO  TO  S00 

CALI.  SPLINE  (SVAL,  YVAL,  DHI,  DDVI , DS, 1 , JTB LS , NDI H , 1 ) 

DO  190  J= 1 , NRC 
X INT=  SH  (J  ♦ 1) 

CALL  SPLINE  (S  VAL  , YV  AL  ,DHI , DDH  I,  DS,  1 , JTBLS,  N DIM  , 1) 

190  HFLP(J) =FTNT 
DO  200  J=  1 , NRC 
20R  YVAL(J)=HELP(J) 
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RETU  RN 

500  CONTINUE 

CALL  SPLINE (S VAL , TV AL ,DW I, DON  I, DS, 1 , JTBLU, * DIB , 1) 

DO  600  J = 1,  NMNLC 
XINT=SWU(J) 

CALL  SPLINE(SVAL,YVAL,DWI,DDWI,DS, 1 , JTBLU, NDIH , 1) 

600  HELP(J)=FINT 

DO  800  J=NLC,NHl 
800  YVAL  (J)=HELP(N-J) 

RETURN 

END 

SUBROUTINE  G JR (A, N, EPS ,N ARON , N ACOL) 

DOUBLE  PRECISION  SOLUTION  OP  A*X=B.  THE  VECTOR  B IS  AUGHENTED  ONTO  THE 

LAST  COLUMN  OP  A.  THE  ANSWER,  X,  IS  ALSO  RETURNED  IN 
THE  LAST  COL  OF  A. 

IPIVOT  IS  A VECTOR  CONTAINING  THE  PIVOTING  ORDER. 

R EAL*8  A (NAROW,NACOL) , D(90)  , B (90) , X (90) , ROWHAX , COLMAX, AWIKOV , RATIO 
$, SUH , DABS , UMAX  1 
INTEGER  IPIVOT  (90) 

I PLAG=1 
NP1=N ♦ 1 

INITIALIZE  IPIVOT, D,B 
DO  10  1=1, N 
IPIVOT  (I)  =T 
ROWH AX=0.0 
B (I)  = A (I,NP1) 

DO  9 J=1, N 

ROWHA  X=0HAX1 (ROWH AX, DABS (A (I, J) ) ) 

9 CONTINUE 

IP  (ROWHAX. EQ. 0.0) GO  TO  999 
D (I)  = R OW H A X 

10  CONTINUE 

GAUSS  ELIMINATION  WITH  SCALED  PARTIAL  PIVOTING 

NH 1=N -1 

IP(NHI.EQ.O) RETURN 
DO  20  K= 1 , NH 1 
J=K 

KP1=K* 1 
I P=TP IVOT (K) 

COLM  AX=D A BS (A  (IP ,K) ) /D (IP) 

DO  11  I=KP1,N 
TP=IPI VOT  (I) 

AWIKOV=DABS(A  (IP,K)  )/D(IP) 

TP(AWTKOV.LE. COLMAX) GO  TO  11 
COLMAX=AWTKOV 
J = I 

11  CONTINUE 

IP (COL MAX . EQ. 0.0)  GO  TO  999 


T FK=I  PI  VOT(J) 

IPIVOT  (J)  *1  PI  VOT  (K) 

IPIVOT (K) = IPK 
DO  20  T*K  P 1 , N 
IP=I  PIVOT  (T) 

A (TP,  K)  =A (IP, K)  / A (IPK  , K) 

R ATIO*-A ( T P, K) 

DO  20  J=K  P 1, N 

A (IP,  .1)  = R ATIO*A(IPK,.l)  *A  (IP,  J) 
20  CONTINUE 
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TF (A ( IP, M) . P0.0. 0) GO  TO  999 

GO  TO  2S 

SET  IFLAG=2  TO  INDICATE  TNABTLITT  TO  FACTORIZE  MATRIX. 

999  IFLAG=2 

H RTTE  (6,  9999) 

STOP 

2S  TF(N.  GT.  1)  GO  TO  90 
X(1)=B(1)/A(1,1) 

RETURN 

90  T P=I PI VOT  (1 ) 

X (1)  =B(IP) 

DO  SO  K=2 , N 
I P=IPTVOT  (K) 

KM1=K-1 

sum=o.p 

PO  UP  J=1 ,RM1 
40  SUM=A  (TP,  J)  *X  (J)  +SUM 
FO  x (K)  =B  (TP)  -sun 


X (N)  =X  (N)  / A (I  P ,N) 

K = N 

DO  70  NP 1MK=  2 , N 
KP1=K 
K=K-1 

IP=I PI VOT  (K) 

SUM=0.0 
PO  60  J = K P 1 , N 
60  SHM=A  (TP,  0)  *X  (J)  *SUM 
PO  X (K)  = (X  (K) -S0K)/A  (IP,K) 

PLACE  ANSWER  VECTOR  IN  LAST  COl  OF  A. 

PO  90  .1  = 1 ,N 
90  A (.7,  NP1)  =X  (.7) 

RETURN 

990  PORN  AT ( 1 HO,  ' ****  ON  ABLE  TO  COMPLETE  L-0  DECOH POSITION  OP  MATRIX') 

END 

SUBROUTINE  DERPS (X , V A LS , RATES) 

RETURNS  DDEL  DX, DUB  DX , DO  TDX  TC  CALLING  ROUTINE.  THIS  IS 

STORED  IN  VECTOR  RATES. 

TBL  COMPUTATION  WITH  PRESSURE  SPECIFIED. 

REAL  K AP,  VALS  (9)  , RATES  (3) 

REAL  A (9  , 9)  , B ( 3)  ,W(9,9),D(3) 

INTEGER  IPIVOT(3)  ,IFLAG 
COMMON/DFR1/DDEL DX, DUBDX,DOTDX,DUIDX1 

COM  MO  N/PE  R2/DEI.T,  UB,0T,UI1  , VT , VB  ,0  DO  I ,T  AOH,  H,  THET  A,  WLST,CFD2, 

$V  ISCOS , N BL 

COMMON/ODE1S/JSTRTS, JENDS , NDIM , S W ( 90)  ,VI  (90) , DVI  (90)  ,DDVI(90) , 

SOS  (90) 

COMMON/ODE 1U/JSTRTU ,JENPU,SWU  (90)  , VI 0 (90) 

COHMON/ODE2UA)TBLU,STBLU(90)  , DSTA  RU  (90),UI1D0  (90)  ,DELT0(90)  , 

$01200 (90) 

rOMMON/SPLTN/XX,UI,DOIDX,ODOI,ISETOP,KMID 
COMNON/TEMPl/XC,IWALLV 
COHMON/TENP2/IEXTT, VFL  (90) 

DATA  KAP/.41/ 

C SET  OP  COEFFICIENTS  FOP  THE  A MATRIX,  ANS  SOLVE  A*BATES«B 

XX=X 

TF (TWALLV . EQ. 1 )GO  TO  100 

CALL  SPLINE  (S  B , VEI.  , DV I ,DDV  I,  DS,  JST  RTS,  J ENDS,  N BIN  , 4): 
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GO  TO  150 

100  CALL  SPLINE (SW0,7IU,DYI,DD7I,DS,JSTRTU,JEWD0,NDIM,*) 

150  C OUT I BOB 

DELT=7  ALS  (1) 

UB*7ALS(2) 

UT=7  ALS ( 3) 

CALL  BL7ALU 
DKU=DELT/  (KAP*UI) 

CALL  TAUHAX  (TAUMEQ) 

TAUH=TAUH  EQ 
A (1,  1)  =TH ETA/ CELT 

A (1,2)= (DELT/OI) *(0.5-0. 75 *7B- 1 . 58949*7T) 

A(1,3) = DK  0* (1 .0-4.0*YT-1. 58949*7 B) 

A (2,1)=0T**2/(KAP*DELT) 

A (2,2) = UT 

A (2,3)  =HT/KAP*UT-U3 
A ( 3,  1)  =1 . O-DELST/DBLT 
A (3,2) =-0.5*DELT/0I 
A (3,3) =-DKU 

B ( 1)  = (K  AP*7T)  **2-2 . 0*DELST*DOIDX/OI ♦THETA/ (XC-X) 

B (2)  =UT*DUIDX 

B (3)  =10.0*TAOH/OI**2-(DELT-DELST*DELT* (7T*0 . 5*  7B) ) *DOIDX/OI 
CALL  FACTOR(A,W,IPIVOT,D,3,IFLAG) 

CALL  SU8ST(W, B,P ATES, IPIVOT, 3) 

RETURN 

END 

SUBROUTINE  PERSEG 

SET  up  x HE  A MATRIX  OF  COEFFICIENTS  TO  SOL7E  LAPLACE *S  ECN 

IN  2-D  USING  PLEMELJ'S  FORB  OF  THE  CIUCHI  INTEGRAL  FORMULA. 
LINEAR  APPROX  FOR  THE  FUNCTION  BETWEEN  NODE  POINTS. 

COMPLEX  C (91 ) 

REAL  XO < 1 20)  ,70(120) 

R EAL* 8 A (86,87) , BB, DT H AG , DREAL 
COMPLEX  Z ERO , ICMPLX 
COHPLEX*1 6 CS  (180)  ,LAHDA0( 180) 

COMPLEX*  16  ZO,ERP,LERP,DMP1,D(!«1,TEHP 
COMPLEX*1 6 CD LOG 

COMMON/EFCV AL/XC (90)  ,YC(90)  ,AL(90)  ,ALNV(90) 

COMMON/GEOM2/N ,NR ,NL , NU , NM 1 , NLC, NRC 
COHMON/PF  SL1/C, A , XO , TO 
DATA  PT/3.  141593/ 

NUPl  =NU*  1 
NM2=N-2 

ZERO= (0.0, 0.0) 

ICHPLX= (0 .0,1 .0) 

C ♦ ♦♦  EXTEND  C ARRAY  ♦♦♦♦ 

DO  30  J=  1 , N 
CS  (J)  =C  (J) 

30  CS  (J ♦ N)  =C  (J) 

C ♦ ♦♦♦  EACH  PASS  CORRESPONDS  TO  ONE  UNKN  ZO  BOUNDARY  POINT  ♦♦♦♦ 

DO  500  M=  1 ,NU 
ZO=CS  (H) 

JSTART=H* 1 
JEND=  NM1 ♦ M 
M J FN  D=J  EN  D-  1 

C ♦♦♦♦  FOPM  C.EOMFTRT  CO BFFICENTS  ♦♦♦♦ 

LAMDAO  (JSTART)  =Z  ERO 
DO  50  J=.T  START, MJEND 
E RP=  (CS  (J  ♦ 1)  - ZO)  / (CS  (.1)  -ZO) 

LE  PP=CDLOG (ERP) / ( ERP- 1.0) 

LAMDAO  (J)  = LAMDAO  (J)  ♦ER”*LERP 
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SO  L AMDAO  {,7  4 1)  =-LERP 
DMP1=CS  (.IStaRT)  -zo 
DMM1=70-CS  (.TEND) 

T FMP=  (DMP  1-DMM  1)  /?.  0 

LAMDAO  (M)  =CDLOG  (DMP1/DMN1)  -(ICHPLX.  *PT) 
f—  (TPM  P* ( 1 . D/DM PI ♦ 1 . O/DMM 1) ) 

TP  (M.  EO.  NPC.OR.M.  EC.  NLC)  L AMD  A 0 (M)  =CDLOG  ( (ICHPLX  * DM  P 1 ) / ( -DBM  1 ) ) 
F-ICMPLT*  (PI/2.0) - (TEMP*( 1.0/DMP1 41 .0/DHH1)  ) 

LAMDAO  (JSTART) =LAMUA9  (.(START)  4 (T  EM  P/D  HP  1) 

L AMD  AO  (.TFNP)  = L AM  DAO  (.TEND)  4 (TENP/DHH1) 

IF  (M. EO. 1 ) GO  TO  70 
M 1 = M-1 

DO  60  .7=1,  Ml 

60  LAMDAO  (J)  =LAM  DAO  (N4J) 

70  CONTINUE 

r 4444  "OEM  A MATRIX  ♦ + »♦ 

u B=0 . 0 

DO  2 60  .1=1,  Nil 
A f M,  J)  =DT  M AG  (LAMDAO  (J)  ) 

DSO  DP=PB  + AL  (J)  *D  REAL  (LAMDAO  (.7)  ) 

A (M,NUP1)  = BB 

04-  AL  (NM  1)  * DRP  AL  (LAMDAO  (NM  1)  ) 

64  A L(N)  *DRE AT.  (LAMDAO  (N)  ) 
sod  CONTINUE 

C 4444  A MATRIX  POP  MU  LATION  COMPLETE  4444 

RETURN 
END 

SUBROUTINE  EFGEOM 

C GIVEN  THE  WALL  LOCATION  COORPIN  AT ES, XW (I) , YW (I) , A LW  (I)  , 

C AND  DTSPL ACFMFNTS  DSTAR  S ( I)  , DSTAR  U (I ),  LOCATE  THE  EOUNDAPY 

C OF  THE  E PC.  NOTE  STBI.S  AND  STBLU  ARE  DISPLACED  ONE  ELEMENT 

r AHFAD  OP  THE  RFST, 

REAL  DSHIFT(OO) 

R F AL*  B A (86,87) 

COMPLEX  C(91) 

REAL  X 0 ( 1 2 0)  , Y0(120) 

COMMON/EFCVAL/X (DC)  ,Y  (90)  , AL  (90)  ,ALNV (90) 

COMMON /G POM  2/N ,NR , NL, NU, NMl, NLC, NRC 

C0MM0N/0DE1S/JSTRTS,  J ENDS  , NOTH  , S W ( 90)  ,NI  (90)  ,DWI  (90)  ,DDWI(90)  , 
SDS  (QO) 

COMMON/O DP1U/.TSTRTU,  JENDU,SWU  (90)  , WI  U (90) 

C 0MM0N/0DE2S/ JTBLS ,STBLS(90)  , DSTARS  (90) , Oil DS  (90) ,DELTS(90)  , 
SUI2DS  (00) 

C0HM0N/0DE2U/JTBLU, STBLU  (90)  ,DSTAPU  (9  0)  ,DI1DU  (90)  ,DEITU(90)  , 
SUT2DII  (00) 

COMMON/PPST.  1/C,  A,XO,  TO 

C0HN0N/SPLYV/XINT,PINT  , DDSTDX, PPPTNT, I SETUP, KM  ID 
C0NM0N/TEHP1/XCMX, IWALLV 
C OMMON/W ALVAL/XW (90) ,TW(90) ,ALW (90) 

NRCP1 =NRC ♦ 1 
N MNL  C=N-N  LC 

IMfTWALLV.EQ.  1)G0  TO  10S 

ON  LOWER  WALL,  DSHIFT (J4 1)  = DSTARS (J) , I.  E.  DISPLACED  ONE 

FLEM  ENT. 

DSHTFT(I)  =D STARS  (N) 

DO  60  J=  2 , NRCP  1 
SO  DSHTPT  (.1)  =DSTARS(J-1) 

C SET  UP  DSTARS  AND  ITS  DERIVATIVE  DDSTS  ON  LOWER  BOUNDARY. 

TINT=SW{1) 

I S PTU  P*0 
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K NID  = 2 

CALL  SPLINE (SW , DSHI FT , DWI , DDW I, DSf 1 , N BCP1 , H DIN , 2) 

X (HI  =XW(N) 

r (N)  =TW(N)  ♦ DSTARS  (If) 

T 0 (N)  =XW(N) 

Y 0 (N)  = YW  (If)  ♦DELTS  (N) 

AL  (N)  = ALB  (N)  ♦ ATAN  ( DDSTDX ) 

DO  100  J=1,NRC 
X INT=  SW ( J+ 1 ) 

CALL  SPLINB(SW,DSHIFT,DWI,DLWI,DS,  1 , N BCP1  , If  DI fl  , 2) 

STNALJ=SI N (ALB  (J)  ) 

COSALJ=COS  (ALB  (0)  ) 

X (J)  = XH  (J ) -DSTARS  (J)  *S  TNALJ 

Y (J)  = YW  (J)  ♦DSTARS (J)  *COSALJ 
X 0(J) =XW (J) -DELTS  (J) *SIN  AL J 
YO  (J)  = YB  (J)  *DE  LTS  (J)  *COS  ALJ 

100  AL  (J)  =ALW  (J) +ATAN  (DDSTDX) 

GO  TO  160 
r 

SET  OP  DSTARO  AND  ITS  DERIV  DDSTS  ON  UPPER  BOUNDARY. 

FLIP  INDICES  TO  HAKE  DSTARO  INCREASE  IN  SANE  DIRECTION  AS  SBU. 

106  DO  110  J = NLC,  NNl 
110  D SHI  FT  (N-J) = DSTARO (J) 

XINT  = SBO  (1) 

TSFTO  P=0 
KNID=2 

CALL  S PI.  IN  F (S  BO,  DSHIFT,DBI,DDBT,DS,1,  NNNLC  ,NDTH,2) 

DO  160  0=1 , NNNLC 
X TNT=SWU (J) 

CALL  SPLINE  (S BU, DSHIFT , DWI ,DD» I, DS, 1, NNNLC, NDIN , 2) 

NNO=N -J 

S INA  L J=SI N (ALW(NNO) ) 

C0SAL0=COS  (ALB  (NNJ)  ) 

X (NNJ) =XW  (NNJ) ♦DSTARO (NNJ)  *SINALJ 
Y (NNJ) = YH (NNJ) -DSTARO (NNJ) *COSALJ 
XO  (NNJ)  =XK  (NNJ)  ♦DELTU(NNJ)  *SINALJ 
Y''  (NN  J)  =YW  (NNJ)  -DFLTU  (NNJ)  * COS  ALJ 
160  AT.  (NNJ)  =ALW(NNJ)  -ATA V (DDSTDX) 

160  CONTI  NOF 
C 160  WRITE  (6,  ROY) 

C RO-7  EOPNATC  1WALI  COORDINATES  AND  EFFECTIVE  FLOW  CHANNEL  LOCATION*//) 
r WP IT  E (6 , B9R) 

C BOR  Porn  AT  (•  0*  ,T4,*  J'  ,T10,  * XW  (J)  • ,T25.*YW(J)  • ,TU0,'  ALB  (J)  • ,T5S, 

C X*X(J)  *,T70,* Y(J)  • ,T86,»  AL  (J)  ' ,T 100 , • DELSTAP'/) 

C WRITE  (6,  900)  (J,XW(J)  , YW  (J>  , ALW  ( J)  ,X  (J)  , Y (J)  ,AL  (J)  .DSTARS  (J)  , 

C $ J = 1 , NPC) 

C WRITE  (6,  900)  (J,XW(J)  , YH  (J)  , A LB  { J)  , X ( J)  , Y (J)  , AL  (J)  , DST  ARU  (J)  , 

(*  X J = NRC”1  , NNl) 

C WRITE  (6, 900)  N ,XW  (N)  , YW(N)  ,ALW  (N)  ,X  (N)  ,Y  (N)  ,AL  (N)  , DSTARS  (N) 

r quo  FnpNATdS,  1P7P16.6) 

RPTORN 

END 


SOPROOTINP  DERSI(T,VALS, PATES) 

C RFTORNS  DDFLDX,DI)BDX,DUTDX,DUIDX  TO  CALLING  ROUTINE,  VIA  THE 

r VECTOR  RATE’S. 

C ^RL  CONFUTATION  WITH  SINOLTANEOUS  ITERATION  RETWEEN  THE 

n BOUNDARY  LAYER  AND  ON  E-  DINFNSIONA  L CORE. 

PEAL  KAO, VALS  (4)  , RA^ES  (4) 

R R AL  A (4,  4)  ,R  (4)  ,W  (4,4)  ,D  (4) 

TNTEGRR  IRTVOT(4)  ,TFLAG 
C ON  NON/OF  R 1 / D D FI  PY  , OUP  PX  , DIITD  X , DUX  D X 


> 


''OMMOM/OHR2/D UR, 'IT, ITT  , VT  , VP  ,fl  PH  I ,TA  UN,  !!,  THETA, TF!  ST,CPU2, 
*VISmS,NBL 

connon/gfon i/xn,  ui,i*h,  udth,x i ,s  i ,sr  nth.costh,  twothp.twothi  , 

ASTN^H  1 ,COS  HI  ,AS,WH,XC,X2,XNAX,XDF 

comon/oofi s/.ts-rts,  j^nds,  ndim,  sh  (9  0)  , hi  (90)  , dhi  (9  0)  ,ddhi  (9o> , 

Sr)  S (HP  ) 

CnHMON/n(lElH/.THT“T'a,J  FHPO,  sun  (90)  ,Hin  (HP  ) 

COMNON/O  DF2S/.7TPLS  , SO’BI.S  (9  0)  , OSTAPS  (9  0)  ,111  1 OS  (9  0)  , OF  ITS  (90)  , 
*91209  (HP) 

C0«»0N/nDF2n/JTBI.n,STBI,U  (90)  ,OSTAPU  (90)  ,(IT10U  (90)  ,DELTU(90)  , 

*11121)1!  (90) 

CON»ON/SPLYN/XX,WTDTH,DHDX  ,ddvpx,isetup,kntd 
C ON*  ON/TF  N P 1/XCN  X , IUALLV 
OATA  KAP/o.41/ 


C HIDT  H=W I DTP  OF  1-D  COPE  SECTION (CHANNEL  SITE  NINOS  BLOCKAGE) 

c gpT  CO^F '■TCI  ENTS  FOP  THE  A NATPIX,  AND  SOLVE  A*PATES=B 


YX  = X 

WO'’’  H R=0 , 0 
CnsTH=1 .0 

IP  (NHT..',0.  1)00  TO  nr 
TF  (YX.LS.  XI  .0P.XX.GE.XDE)00  TC  IOC 
CO  S'r’H  =COS  TH 1 
T WOTH  R =T  HO  T H 1 
I""  CONTTNHE 

PSLT=VALS  (1) 

HB=VATS(2) 
tJT  = V M.  S ( 1) 

H I - V A L S (9  ) 

U toth=ot 
0 HD Y = 00ID  X 

TF ( (HB-OT) *IJT  . LF. 0,0)  GO  TO  110 
*77’  = — YTT 

1 ^ r ALL  PLVALT! 

nifn=nF!  t/  (kap*ijI) 

CALI.  TAON  AX  (T  AUN  RQ) 

TAHN  = TAUNFO 

■~ THALLV  = r IS  TH"  lOUEP  H ALL . IH  AL  IV=  1 IS  UPPER  HALL. 

TF  (TWALLV.FO. 1 ) GO  TO  130 

CALL  SPLINE  (SH, HI  , DWI , DOW  I , DS , JSTRTS, JENDS, ND IH , 4) 

GO  TO  140 

130  C4LL  SPLINE  (SHH,HIfI,  DHI,  DDW I , DS,  JST  PTO,  .1 ENDW,  N DIN  , 4) 

140  CONYTNU'' 

ODSTDX  = 0ELST*DDPLDX/DFLT-DELT*  (VT* ,5*VB)  *DniDX/UI*DEIT*DUTDX/ 
1 (KAP*OT)  ♦DELT*D"BDY/(2.o*0T) 

TF  (ABS(PDSTDX)  . GT  . 1 .E-4)  XCHX=  (0  . 5*  H H- DELST)  /DDSTDX 
A (1,  1 ) =THETA/DELT 

A ( 1,  2)  = (DELT/HT)  * (.S-0.75*  VB-  1. 58949*VT) 

A (1,  3)  =DKt1* (1 .O-4.0*VT-1 ,S89  49*VB) 

A (1,4)  =2 . 0*nEL ST/T1 1 
A (2,1 ) =0T**2/  (KA  P*DELT) 

A (2, 2) =DT 

A (2,3) =UT/KAP»MI-HB 
A (2,  4)  =-I1T 

A (3,1)  *1.  O-peIST/DELT 
A (3,2)  =>0  ,«1*D®LT/I1T 
A (3,3) =-PK0 

A (3,  4)  = (DELT- DELST  *DELT*  (VT*0.S*  VB)  )/W 
A (4,1) =A  (3,1)  -1.0 
A (4,2) =A (3,2) 

A («,  3)  =-DKIJ 

A (0,4)  ='(COSTH/(NBL*OI))*(HIDTH-NBL*DELST/COSTH)  ♦ ( D ILT/O  2)  • 


U45 


* (OT/!UP*O.S*UB) 

R (1)  * (K  AP*TT)**2»THRTA/XCHX 

B ( 2)  *0.0 

B (3)  * 1 0.0*TAUR/UT**2 

B (4)  *-DWDX*COSTH/NBL 


C 

r 

C IF  OT  BECOBES  SB  ALL,  THEN  THF  BATBIX  BFCOBES  ILL-CONDITIONED. 

C FREEZE  DDTDX,  AND  REBOVE  THE  DSF  ECN.  SOLVE  THE  REDUCED  SET. 


IF  (ABS  (TIT)  ,GT. 0.025)  GO  TO  200 
A (1,3)=A(1,4) 

A (2,  1)  = A(3,  1) 

A (2,2)  =A  (1,2) 

A(2,3)=A(3,4) 

A (3,1)  *A  (4,1) 

A ( "),  2)  *A  (4,  2) 

A (3,3)  =A  (4,4) 

B (2)  = B (3) 

B (3)  ^ b (4) 

CALL  FACTOR  (A,  W,  IPIVOT.D,  3,IELAO) 

CALL  SUBST(W,B,RATES,IPIVOT,3) 

RATES  (4)  = RATES  ( 3) 

RATES  (3)  =DUTDX 
RETURN 
C 
C 

200  CALL  FACTOR  (A, W,IPIVOT,D,4, IF LAG) 

CALL  SUBST  (W,  B, PATES,  IPIVOT.4) 

RETURN 

END 

SUBROUTINE  SPLINE (X,P,FP,PPP, DX, JST ART,  J END,  NOTH, TNTERP) 
C 

C ♦♦♦♦  SPLINE  IN  TENSION  FIT  OF  P(X) 

C ® IRST  DERIVATIVE  AT  POINT=FP 

C SECOND  DERIVATIVE  AT  POTNT=FPF 

C START  AND  END  OP  INTERVAL=.1  START,  J END 

C TENSION  FACTOP=SIGB  A 

C ROUTINE  PT  R TNEHART 

c 

REAL  X(NDTB)  , F(NDTB)  ,FP(NDIB)  ,PPP(NDIB),DX(NDIH) 

REAL  A (101)  ,B  (101)  ,C(101)  ,D(101) 

COBBON/SPL YN/XINT,FTNT,FPINT, FPPINT,TSETUP,KBID 


C 

SOL  I NE  FTT  np  F (X)  USED  POP  FINDING  FIRST  S 2ND  DERIVATIVES  AT 

C THE  POINTS,  r,  ALSO  FOP  TN""FRPOT  ATTON.  ISET0P=  CCUKTER  CF  # OP 

C TIRES  ROTTTINE  CALLED  USING  SAHE  PIT.  WHEN  ISETUP=0  FPP  TAEIE 

C DETFRNINED.  CUPTC  PtlNOUT  END  CONDITION.  KB  T D IS  GUESS 

C INDEX  OP  INTERVAL  WHEPE  X(FPTD)  < X INT  < X(RRTD*1).  WHEN 

C (INTFPP=1,  FIND  EINT),(=2,  FIND  EPTNT),(=3,  FIND  FFPINT)  , 

C (=4,  FIND  PINT  e EPINT) . EOR  INTPRP>4  NO  INTERPOLATION,  FIND 

C ONLY  DERIVATIVES  AT  POINTS.  FOR  INTERP=0  SETUP  ONLY. 


STONA*2.  S 
SS=SIORA*SIGRA 

S IGB  A =STO  R A*  (.1  EN  D-.J  START)  / (X  ( JEND)  - X ( JSTA  RT)  ) 

STSO^I.O 

I®  (TSF'"UP.NE.0)  GO  to  1 BO 
JFNO  B 1 =,1E  ND  -1 
DO  120  ,1  = .T START,. 1®NDB1 
1?0  DX  (.1)  = X GUI)  -X(.l) 

DX  (.TEND)  - OX  (.1  END-  1) 

,1ST"1  =.Tr’TART*  1 


> 


I 


no  mo  i=jstp i,.ifnohi 

H 1=DX  (.1) 

HiiinxiJ-') 

SH1  = STONA*H1 
*=: Hf!  1 = SIGH  A*HH  1 

A (.T)  = ( (1.  "/HH  1)  - (SIGH  A/S  INH(SHHI))  ) *SISQ 
BP1=SIGNA*  ((1  ."/TANH(SHHl)  )♦(■».  O/^ANH  (SHI)!) 

p (j)  = (BPi-  <1.  o/p  i)  - <1  .o/hhi)  ) *siso 

C (J)  = ( (1  . O/HI)  - ( SIGH  A/SI NH< SHI)  ) )*STSQ 
0 (.1)  = ({”  (.1*H  -F(J)  )/H1)-((P(J)  -F(J-I)  )/HHl) 

14"  CONT  INU  p 
JOPPS  = " 

NnTAG=101 
A (.7START)  = ".  " 
r (JFND) =0 

CUBIC  OPN011T  END  CONDITIONS 

.1  =JST  APT 
HT=DX  (J) 

SHT  = SIGN  A*DX  U) 

B (J)  = ( (SIGMA*2."/TANH  (SHI)  ) - (2.0/HI))  *S  IS0»B  C J ♦ 1 ) 
C (J)  = (1  . "/HI- SIGH  A/SINH(  SHI)  ) *SISQ+C(J*1) 

>1  = JFND 

Hl=nx  (J) 

NHI=SIGNA*DX (J) 

HMT="T  (.T - 1 ) 

A (J)  = (1.0/HM1-STGHA/STNH  (HMT)  ) *SISQ*A  (J-1) 

" (.1)  = ( (STGMA*2.P/TANH(SHT)  )-(  2.0/HI))  *S  ISCK  B (J- 1) 
C (JFND) =0.0 

CALL  TPT0  AO  (JSTP1  , JEN  OBI  , A , R,  C,  D , FPP,  NDIAG.N  DIH) 
p PP(.T  ST  APT)  =0.0 
PPP  f JE*»P)  =0.0 

IF  (TNTFRP.LF.O)  00  TO  1000 
IP"  IF  (INTEPP.GT. 4)  GO  TO  700 


IF (TNTFRP.LF.O)  GO  TO  1000 
IF (INTFPP.GT.  4)  GO  TO  700 

--FTN0  TNTBPVAL  OF  INTERPOLATION.  X(KHID)  < TINT  < X(KtlID< 
IF  (X  ( JS^ART)  . GT.  X (.TEND))  GO  TO  250 

TF(XINT.  LE.  X (JSTAPT)  . OP.  X I NT . GE.  X ( J END)  ) GO  TO  2000 
- * X IS  A NONOTONICALLY  INCREASING  PUNCTTON  BITH  THE  INDEX . 


X(KNIDFl) 


200  IF  (XINT. GF.X (E H ID) ) GO  TO  220 
ENID  = EHTD-1 
GO  TO  200 

22"  IF  (XTNT.LF.X  (ENTD*1) ) GO  TO  300 
KNTD  = ENID* 1 
GO  TO  220 

y TS  A NONOTONICALLY  DECREASING  FUNCTION  BITH  THE  INDEX. 

250  IFfXINT.GE.X  (J START)  .OR.XINT.LE.X(JEND))  GOTO  2000 
25"  TMTTNT.LP.X(KNID)  ) OO  TO  270 
KNTD  = EHID-1 
GO  TO  260 

27"  TF  (XTNT.GE.X(EHIP«-1)  ) GO  TO  300 
EMI"  = ENID* 1 
GO  TO  270 

PERFORM  INTERPOLATION. 

3""  OELX  = XIVT-X(KHID) 

EHTDP1  = EH  IDE  1 
DEI.XP  = X (KHIDP1) -XINT 
"YEN  ID  = OX(EHID) 

GO  TO  (400,500,600,400),  INTERP 

INTERPOLATE  FOR  F(XINT)  = FIBT 

4 " 0 P TNT  =P  PP (EMI D) *STSQ*SINH(SIGHA  *DELX P) /S IN H (SIGH  A* DIEM  ID) 

F.  ♦ ( F(E  HID)  -FPP  (EH  I D)  *S  ISO)  * (DELXP/DX KHID) 

F ♦ (FPP  (ENIDP1)  *SISQ)  • (SINH  (STGHA*DELX)  ) /SINH  (SIGHA*DX  (BID) 
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6*  (P(KHIDPI)  -FPP (K HI  DPI)  *SISQ)  *DELX/DXKHID 
C.O  TO  (1000,  *500,600,  500)  , INTFBP 

C INTERPOLATE  POP  FP(XINT)  = FPINT 

500  FPINT=-SIGHA*FPP (KHID) *COSH  ( S IGB  A*  DELX  P) /SIHH (S IGB A*DXKH ID) 
*-  (F(KBID)  -FF  P (KB  ID))  /DXKHID 

$*SIGHA*FPP(KHIDP1) *COSH (SIGHA*DELX) /S INH ( SIGB  A*  DXKHID) 

*♦ (P(KBIDPI) -FPP(KHIDPI))  /DXKHID 

GO  TO  1000 

C INTERPOLATE  POF  FPP(XINT)  = PPPIBT 

600  FPPTNT=SS  * (FPP  (KB  ID)  *SIHH  (STGHA*DEL  XP) /SIHH (S IGHA*DX RBID) 
*»PPP  (FBI DPI) *SIBH(SIGBA*DEIX) /SI MH ( SIGB A* DXKB I D) ) 

GO  TO  1000 
700  COBTIBOE 

r ♦♦♦♦  INTERPOLATION  FOR  PP  AT  POINTS,  AVERAGE  OP  FORWARD 

C AND  BACKWARD  FOEBOLAS  ♦♦♦♦ 

DO  750  J=JSTP1,JENDH1 
C ♦♦♦♦  BACKWARDS  DIFFERENCE  ♦♦♦♦ 

F PB=P  PP (J -1 ) * (~SIGHA/STNH(SIGHA*DX (J-1))) 

B-  (F(J-1)  - PPP(J-I)  ) /DX  (J-1) 

B+PPP (J) *STGBA*COSH (SIGBA*DX (J-1) ) /S IN H (SIGB A* DX (J-1) ) 

B*  (F  (J)  -FPP(J)  )/DX(J-1) 
r ♦ ♦♦♦  PORWARD  DIFPFRBNCE  ♦ ♦♦♦ 

FPF=PPP(J) * ( (-SIGBA)  *COSH(SIGHA*DX ( J)  ) /S INH ( SIGB A* DX  (J) ) ) 

B-  (F  (J)  -PPP  (J)  ) /DX  (J) 

B*  P PP  (J*  1)  *SIGBA/  (S  INH  (SIGBA*DX  ( J) ) ) 

B + (F(J*1)-FPP(J+1) )/DX (J) 

C ♦♦♦♦  AVERAGE  FORWARD  AND  BACKWARDS  DIFFERENCES  ♦♦♦♦ 

F P (J)  =0.5*  (FPF* FPB) 

750  CONTINUE 

C OSE  FORWARD  DIFF  FOR  START  SEGBENT,  AND  BACKWARD  DIFP 

C FOR  THE  LAST  SEGBENT. 

J=JSTART 

FP  (J)  =PPP  (J)  * ((-SIGN  A)  *COSH(SIGBA*DX(J))/SINH  (SIGBA*DX(J)  ) ) 
R-  (F  (J)  -PPP(J)  )/DX(J) 
fi+FPP(J*1) *SIGB A/(SINH (SIGHA*DX (J)) ) 

R ♦ (F  (J  ♦ 1)  -PPP(J*1)  ) /DX  (J) 

J =JEN  D 

F P ( J)  =PPP  (J-1) * (-S IGB  A/S INH(SIGBA*DX(J-1) ) ) 

R-  (P(J-1)  -FPP(J-I)  ) /DX  (J-1) 

B**P«  (J)  *SIGBA*COSH  ( SIG BA*DX  ( J- 1 ) J/SINH  (SIGH  A*  DX  (J-  1)  ) 

&♦  ( F ( J)  -PPP  (J)  ) /DX  (J-  1) 

IF  (TNTERP. LT. 5)GO  TO  1000 
DO  BOO  J=J START, JEND 
BOO  FPP(J)  =FP P (J)  *SS 
1000  ISFTUP  = I SETUP* 1 
RETURN 

2000  JOPPS  = JEND 

IF(ABS  (XINT-X  (JST  ART)  ) .LT.ABS  (XINT-X (JEND)  ) ) JOPPS  = JST ART 

C USE  PORWARD  DIFF  FOR  START  SEGBENT,  AND  BACKWARD  DIFF 

C FOR  THF  LAST  SEGBENT. 

J = JST  APT 

»P(J)  =P“P  (J)  * ( (-SIGN A)  *COSH  (SIGBA*DX  (J) ) /SINK  (SIGBA*DX(J)  ) ) 
B- (F  (J)  -FPP(J)  )/DX(J) 

B+FPP ( J*1 ) *SIGHA/(STNH(SIGBA*DX(J))  ) 

6*  (P  (J  ♦ 1)  - FPP  ( J ♦ 1)  )/DX  (J) 

J =JEND 

FP(J)  =FPP  (J-1) * (-S IGBA/STNH(SIGBA*DX(J-1)  ) ) 

B-  (F  (j-1)  - FPP  (J-1)  ) /OX  (J-1) 

B ♦ FPP  (J)  *STGNA*COSH(SIGBA*DX(J-1)  ) /SIHH  (SIGH  A*  DX  (J-  1)  ) 

6*  ( F ( J ) -FPP  (J)  ) /DX  (J-  1) 

PINT  = P (JOPPS) 


FPTN'’’  = ”P(JOPPS) 

FPPTNT  = FPP(JDOPS) 

R ETU  ° N 
END 

pnNCPION  YINT  (X, Y,XINT) 

r GIVEN  3 COORDINATES  (X,Y),  FIT  SECOND  ORDER  LAGRANGE 

C POLYNOMIAL  AND  RETURN  THE  VALUE  YINT,  CORRESPONDING  TO  TINT. 

RPAT,  X(3)  , Y(3)  ,X  INT 
P 1=X  (2)  -X  (1) 

P2=X  (3)-X  (2) 

D3  = X (1)  -X  (3) 

YTNT  = -Y  (1) * (XINT-X (2) ) * ( X T NT-X ( 3) ) / (D 1 *D3 ) 

* -Y  (2)  * (XINT-X  (1)  ) * (XTNT-X  (3))  /(D1*D2) 

♦ -Y  (3)  * (XTNT-X  (1)  ) * (XINT-X  (2) ) / (D3*D2) 

RETURN 

END 

SUBROUTINE  STAND 

C GENERATE  NODE  POINTS  FOR  STANDARD  DIFFUSERS.  STRAIGHT  WALLED 

C UNITS  WITH  BOTH  WALLS  D IV  ERG  I NG  (GEOB  T=  • STDD' ) OR  ASSYHHETRIC 

C UNTTS  WT^H  ONE  DIVERGING  WALL(GEONT= ’HALF') . FOR  NOHENCLATURE 

C SEE  USERS  GUIDE. 

C 

C 

REAL  N,L 

COMMO N/BLT V/H S ,DELSTS , HU ,DELSTU 
COM*ON/DER1/DDFLDX,DUBDX,DUTDX,DUIDX 

COMMO  N/DE  P2/DELT, UB,UT,UI  , VT , V B ,U DUI ,TAU H, H, THETA, DELST,CFD2, 
SVTSCOS.NBL 

COMBO N/EFCVAL/X  (90)  , Y (90)  , AL  (90)  ,ALNV  (90) 

COMMON /GFOM1 /N,W1 ,TH, WIDTH, XI , B1 , S INTH, COSTH, TWOTHR ,TWOTH1 , 

IS  TNTH1,COSTH1,AS,WH,XC1,X2,XHAX,XDE 
C 0M**O  N/G  EOM2/NS,NP,NL,NU  , NSB1  , NLC,  N RC 

COMMON/ODE  IS/ JSTRTS ENDS,NDIM,SW(90)  ,WI  (90) ,DWI  (90)  ,DDWI(90>  , 

$DS  (90) 

COMMON/O  DEI U/JSTRTU , J ENDU , SWU (90) ,WIU (90) 

COMMON/IN  IT  AI./DELST  1 , H 1 , U1 1,  IPR  1 
C OM MO  N/N  STD/I  D1  , I D2  , I 03  , NST,  SWT  (90  ) 

COMMON/PR INT/ IPP, NORM  PR, CPERO R, ITH  AX 
COMM ON/S PL YN/XX, WT,DWT,DDWT, ISETUP, KHID 
COMBON/TEMP1/XC,  IK  ALT.  V 
COMMON/WALVAL/XW (90) , YW(90) ,ALW (90) 

DATA  PT/3.  141593/ 

READ  (5,402)  XI  ,RC1 , N ,RC2 , X2 , W 1 , TW OT HD,  AS 

C IF  BOTH  THF  INLET  AND  DIFFUSING  SECTIONS  ARE  OF  ZERO  LEIGTH,  QUIT. 

TF(N. EO.0.0.AND. X1.FQ.0.C) RETURN 
WRITE  (5,903) 

WRITE  (5,904)  X 1 ,RC  1 , N , RC2 , X 2 
WRITE  (5,910) 

W RITE  (5,  0 15)  W 1 ,TWOTHD,  AS 
r 


--STORE  STARTING  VALUES. 

IF  (AS .LF.o.O)  A S=R .0 

READ(5,901)N1,NC1,N2,NC2,N3,ND1,ND2 
WRITE  (5, 9 34) 

WRTtr  (f,935)N1,NC1,N2,NC2,N3 
R EAD ( 5 , 900)  B 1, HI , VISCOS , XC 
TF  (XC.EO. 0.0) XC=1 . E5 
WRITE  (5,920)  B1,UI , VISCOS, XC 

-- N,W 1 , DFLST  (FT) , TWOTH  (DEGREES) , BI(N-D),  UI (FT/SBC)  , VISCOS ( FT2/SEC) 
READ  (5,905) HS , DELSTS , HU, DELSTO 
RE  AD (5, 405) TPP.NORMPR, ITHAX,CPEROR 
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H = HS 
H1»HS 

DELST-DSLSTS 
DELST1 =DELSTS 
THETA*PELST1/H1 

rr (ho. he. o.O) go  to  so 
Ht1  = H1 

T)ELSTO=DELST1 
60  CONTENDS 

W RTTB  (6, 007) HS,DELSTS , HO ,DELSTU 

WRITS  (6.90R)  IPR,NORHPR#ITHAX,CPEROR 

NRC=N1*NC1*N2*NC2+N3 

NLC=NI»C*1 

N S = 2*  2*N  R C 

W RCP1 = NRC* 1 

NSM1=NS-1 

NSH?=NS-2 

N END=0 

TWOTH  P=TWOTHD*PI/180 .0 
THP=TWOTH  P/2. 0 
T AND2  =TAN  (THP/2.0) 

'r  A N^H  1 =TA  M (THP) 

S INT  H 1 = SIN  (THR) 

COS^HI  =COS  (THR) 

TWOTH  1='’’ WOT  HR 

SINTH=SINTH1 

COSTH=COSTH1 

PC  1 MT=RC1 *TAND2 

PC2MT=RC2*TAND2 

C 1 = X 1-PC  1 HT 

C ?=X 1 ♦ RCl  NT*COSTHl 

C 3=X  1 ♦ N-RC2HT  *COST  H 1 

C4=X1*N*RC2HT 

TT  T 1 = 0 I 

H 1=H 

DBLST1 =PE  LST 
T PR  1= IPP 


C W 1 , W 2 ARE  INLET  , EXIT  WIDTHS (FT) ,L  IS  SLANT  LENGTH  AtONG  WALL. 

' ’ L = N/COS^H1 

X DE=X  1 *1. 

tmax=xoe*x2 

W2=W1 »2. 0*L*SI NTH  1 

C STARTING  VALOE  AT  WODE  ZBRO(=NS). 


X (NS)  =o.o 
T (NS)  =0.0 
A L (NS) =0. 0 
X (NSH 1 ) =0  .0 
T (NSN 1 ) =W1 
Will)  = W1 

sw  (1)  =0.0 
swo (i) =0. o 

AL(NSNl)  =0.0 

r rooPOTNATES  POP  INLET  SECTION,  N1  SEGMENTS. 

X 0=T  (NS) 

T F ( N 1 .PO. 0) GO  TO  120 

c APT'T,H=ETIC  PROGRESSION  FROM  INLFT  TO  THROAT. 

IF(NDl.FO.r')NDl  = S*Nl 
Pt=Tl-RC10T 
A = PL  /N  1 

0*0.  f' 
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n o TO  T TO 


t 

i 

4 


IE(N1 . FO. 1)00  TO  100 
a =fl/ndi 

0 = 2.  0*  ea-A*N  1)  / (N  1*  (N  1-  1)  ) 

NSTAPT=o 
N ENP  = NS‘,,ART*N  1 

loo  no  no  .1  = 1, Ml 

0X  = A*  (N1- J) *0 
X (J)  =X0*DX 

Y (J)  =0.0 

» Lf.1)  =0.0 

SW  (J*1)  =X  (J) 

i io  xo=x  (j) 

n 

r THROAT  CURVE,  NCI  SEGMENTS. 

1?o  r=  («(C1  .E0.'>)  GO  TO  140 
XCl=PC1«T*  (1. O+COSTHI) 

OX rXCL/VCI 
N STA  RT=N 1 ♦ 1 
VFNn  = »)STA  FT*  NCI  -1 
00  130  J=N ST ART, NEMO 

x (.i)  = xo*nx 

Y (J)  = -PC  1 + S OPT  (PC1**2-  (X  (J)-  (X1-RC1HT))  **2) 
A l (.1)  =-APSTN  ( (X  (J)  -Cl)  /RC1) 

SH  (.1  + 1)  = SW  (N START)  *RC  1 *A  B S ( A L ( J)  ) 

1 10  XO=X (J) 


DIFFUSING  SECTION,  N2  SEGMENTS. 

14°  IF (N2 . PO. 0) GO  TO  160 

Apt  THMFTTC  PROGRESSION  FROM  THROAT  TO  TAILPIPE. 

IF  (ND2.  p0 . 0)  ND?=S*N2 
FL=N- (PC  1 MT+PC2MT) *COSTH 1 
A=  FT  /N  02 

0=2.0*  (FL-A*N2)/(N2*(N2-1)  ) 

V STA  R T = NE  ND* 1 

N FND=  N S’r  ART*N  2-  1 
PO  ISC  J=  N STA RT, N ENO 
DX  = A ♦ (J-NST  ART)  *0 
X (J)  =XO+PX 

Y (J)  = (Xl-X  (J)  ) *t A N TH 1 
AI.U)  =-THR 

SW  (.1*  1)  =SW  (NSTAPT)  ♦ (X  (.1)  -C2)  /C0STH1 
ISO  X0=X  (.1) 


TAILPIPE  INL  Fi*  CURVE,  NC2  SEGMENTS. 

1«0  TP (NC2.E0. 0)GO  TO  1R0 
YTFMP  = Y (NEND) 

XCI  -RC2NT*  (1  . ONCOST  HI) 

0X=XCL/NC2 
N ST  ART  = N EN  D*  1 
NEND  = NST A RT*  N C2- 1 
DO  170  .1  = NST  A RT , N EN  D 
X (.1)  =X0*PX 

Y (J)  =RC2-N*TANTH1-S0PT  (RC2**2-  (X  (J)  - ( XI  *II*PC2BT)  ) **2) 
D 7 =S  Q RT  f(X{J)  -Cl)  **2*(Y(J)-TTEHP)**2) 

RETA  = 2.0*ARSIN  (DZ/ (2 . 0*RC2) ) 

SW  (.1*1)  = S«  (NSTAPT)  *PC2*BETA 
A L ( .1 ) =PET  A-THR 
no  XO  = v (J) 

TATLPTPF  SECTION,  N3  SEGMENTS. 

1R0  I E (N 3 . FO. 0) GO  TO  200 
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YCL=X2-RC2NT 
DX=XCL/N3 
N STA  RT=NE  ND  + 1 
NEND=NSTART*N  1-1 
DO  190  J=NSTART,NFND 
X(J)=XO»DX  1 

T (J)  =-N*TANTH1 
At  Ml  =0.0 

SW  (J*1)=SR(N  START)  *X  ( J)-C4 
190  X0=X (J) 

HAP  UPPER  BOUNDARY  PROH  LOWER  WALL. 

20D  no  250  J=1,NRC 
NSH1 H J=NS  H 1 -J 
X (NSB1JM)  =X(J) 

Y (NSB1HJ)  =W  1- Y (J) 

AL  (NSB1HJ)  =-AL  (.1) 

SWU(J*1)  =SH(J*1) 

25"  HI  (J4-1)=W1-2.0*Y  (J) 

J STPTS  = 1 
J ST  RT  U = 1 
JENDS  = NPCP1 
J ENDfJ  = NPCP  1 
J END  PI =JENDS* 1 

W RTTE  (6,  1212)  (J,X(J)  ,Y(J),AL  (J)  ,WT  ( J)  ,SW  (J)  ,J  = 1 ,JENDS) 

1212  PORBAT  ( ' ’,15,1 P5E15 . 5) 

DO  300  J=1  ,NS 
XV (J) =X  (J) 

YW  (.1)  =Y  (J) 

300  ALW  (J)  = AL  (J) 

VPTTE  (5,1  313)  (J, X(J) ,Y(J) ,AL( J) , J=JENDP1,NS) 

1313  POPBATC  * , 1 5 , 1P3P15.5) 

r 

C IP  ID1  = U THEN  SET  HP  STANDARD  DTFPOSER  WITH  ONLY  1 DIVERGING 

C WALL,  AT  AN  ANGLE  -T HET A=- ( TWOT HD/2) . NOTE  TVOTHD  ENTERED  BUST 

C BE  DOUBLE  THIS  VALUE.  HOPTPIES  OUTPUT  FROH  STAND  BY  CHOPPING 

C OFF  TOP  WALL. 

I F (TD 1 . NE . 4) RETURN 
NSTAR'"=N1  »2 
SST=SW  (N  1 ♦ 1) 

DO  350  .1=  NSTA  RT,  N RCP1 
SWU(J)=X(J-1) 

VT  (J)  =WT  (J)  - (BI(J)-WI)  /?  . 0 
N 1 J = NS-J 
X W (NBJ)  =XW  (J-  1) 

YW  (NB.1)  =W1 
3 50  AI.W  (NBJ)  =0.0 
DO  400  J=1 ,NS 
X M)  - XW  (J) 

Y (.1)  = YW  M) 

400  AIM)  = ALW  (J) 

THD=TWOTHD/2. 0 
WPITE  (*,049) THD 

W RTTE  (5,  1515)  (J,x  M)  ,Y(J)  ,AL  ( J)  ,SW  (J)  ,.1  = 1 ,JENDS) 

1515  POPBATC  ,,T5,1P4F15.5) 

WRTTE  (5,1414)  (J  , X (.1)  ,Y(J),AL(J)  , J = .1ENDP1 , NS) 

1414  POPBATC  ',T5,  1P3F15.5) 

RFTURN 

940  FORB AT »• 1***STANDARD  DTPPUSER  WITH  1 DIVERGING  WALt  AT  AN  ANGLE • 
9,p7.1, * (DEGREES)  ****//•  WAIL  COO  ED  I NA TES-NO DE# • ,T30, • X-COORD* , 

5 T 4 5 , • y-roORD*  ' , ' ALPHA  (PAD)  • ,T7r>,  'W  IOTH'/) 
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<100  PO  PM  AT  (4  F 1 0 . 5 ) 
o,->1  POEMA"  (ij  10) 
dm  nOPMAT  (RF 1 0. 5) 

grp  FORMAT  (' -DTFFUSFO  GEOMETRY-INLET  ( X 1 - FT)  , THROAT  RAH  (RC  1-  FT)  , DIPF 
♦USING  LENGTH  (N-FT)  , EXIT  P AD  I OS  (PC  2- FT)  , T AI  I.  PI  PE  ( X2 -FT)  • ) 

DO  4 pop*  a T ( ' ' ,T1B,P1o.5,T34,P10  . 5 ,T54 , F10,  5, T78,  pi  o.  5, T99,  P10. 5//) 

OOP  PQRMAT  (UFIO.O) 

OPS  FORM  AT(DT10,E1O.O) 

RDR  format  <• 0 INLET  PL  V A LO FS : LOW KR  W ALL -H3 « , F5 . 2,  • , DELSTS3 • , E 1 2 . 5, 

%’  (FT)  */T 17,* UPPER  WALL-H3' ,P5 .2,  *,  DELSTU3  * , E 1 2 . 5, • (FT) • //) 

DOR  PORMATC-RL  PRINT  INTE  RV  A L ( I PP)  = • , I 2,  • , NORMPR3*, 

STD,',  *AX  # ITERATIONS3* ,12, ' , MAX  ALLOWABLE  CP  ERROR3* , 
♦1°F12.S//) 

Dio  pop* at  (•-* ,T1 f ,' WIDTH (W1-FT)  , TW CTHD { DEGREES)  , ASP FCT-RATIO* ) 

DIR  FORMAT  ' ,T1 B,F10.5,T  14,p10. 5,T54,F1 0.5//) 

D?*  FOPMAT(1HO,'  B1 ,01 , VISCOS, XC= • ,2F12 .5, FI  2. 7, E12. 5/) 

otu  PORM  AT  ('  S PGM  FIT  DISTRIBUTION  - INLET, THROAT  CURVF,  DIFFUSING  SECT 
^T I ON,  PX IT  CURVE,  TAILPIPE') 
ddr  FOPmaTC  ' ,T24,I2,T31  ,T2,T45,  T2,T6  4,I2,T76,  12//) 

FND 

SUBROUTINE  DERSTI.  (X,VALS, RATES) 

C RETURNS  DDELDX, DOBDX, DUTDX,DUIDX  TO  CALLING  ROUTINE,  VIA  THF 

C VPCTOP  PATES. 

C TP  L COMPUTATION  WITH  SIMULTANEOUS  ITERATION  BETWEEN  THE 

C BOUNDARY  LAYER  AND  COPE,  ASSUMING  LINEAR  VELOCITY  VARIATION 

C BETWEEN  THE  0 PnE R AND  LOWER  DEI.TASTAR  LINES. 

REAL  K AP,  VALS  (4)  ,RATES(4) 

REAL  A (4,  4)  , 9 (4)  , W (4,4)  ,D  (4) 

TNTFGRB  I PTVOT  (4)  ,TFI.AG 
COMMON/OER  i/udeldx , dubdx ,dutdx, duidx 

C0MMON/DEP2/DELT,UB,UT,U I ,VT,VB,UDUI,TAUH, H, TH ETA , D ELST,CPD 2, 

♦ V ISCOf.NBL 

COMMDN/GFOM  VXD, W 1,TH,WDTH,X1 ,B1 ,SINTH,COSTH,TWOTHR , TWOTH 1 , 

IS INTH1 ,COSTH1  , AS,WH,XC1,X2,XMAX,XDE 
C0MMON/I.INEAR/WDTE(D0)  ,DU2D  (DO)  ,DDU2D  (90)  , W MD , D WED  X , UEFF , CU EEX 
COMMON /ODF1S/J ST RTS, JEN DS,NDIM,SW(90)  ,W I (90) , DW I (90) ,DDWI (90)  , 

*DS  (90) 

C OMMON/OD E 1 U/JSTRTU, J FND U, SWU  (90)  ,WIU(90) 

COMmon/ODF2s/JTBLS,STBLS (90)  , PST APS (9  0) ,UI10S (90) , BELTS (90) , 

*U I20S  (DO) 

C OMMON/ODE2U/ JT  BLU  ,ST  B IU  (90)  ,DSTAPU(90)  ,UI1OU(90)  ,DELTU(90)  , 

♦DI2DU  (90) 

COMMON /SPLYN/XX,Y,DYDX,DDYDX, ISFTUP,KMID 
C OMM  ON/TE  MP1/XC,TWALLV 
DATA  KA n/0.41/ 


C WT  CONTAINS  (WDTF-DST ARU) . WTU  CONTAINS  UI2DU(WITH  PROPER  INDICES). 

C SET  COFFFICTPNTS  FO»  THE  A MATRIX,  AND  SOLVE  A*RATES=B 


XX=X 

DEI.T=VALS  (1) 

1J  P=V  ALS  ( 2) 

UT=V  A LS O ) 

U T=V  ALS  (4) 

Y =UT 

D YDX=  R AT  ES  (4) 

IF  ( (IJB-UT)  *UT.LE.O  .0  ) GO  TO  110 
UTs-UT 

110  CALL  BLVALU 

DK0  = DELT/  (KAP*UT) 

CAT  L TAUMAX  (TAUMEQ) 

T AUM=TAUMEQ 

C IWALLVpO  is  THF  LOWER  W A Lt. IWA LLV 3 1 IS  UPPER  WALL. 
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CULL  SPLINE  (SB , HI  , DHI  , DDHI  ,DS , JSTRTS, JENDS, NDIM , 4) 
BHT)=Y 

DBEDX=DYDX 

BEPP=HHD-DELST 

CALL  SPLIHE  (S B ,WI  0,  DU2D,  DDU2D, DS, JSTRTS, J 2NDS, NDIM , 4 ) 
UEFF=0.5* (T  *0 I) 

DUEDX=DYDX 

A (1,1) =THETA/DELT 

A (1,2)  = (DELT/UI)  *(.  5-0. 75* VB-1 .5894  9*  VT) 

A (1,  3)  =DKU*  (1  ,0-4. 0*TT-1. 58949 *VB) 

A (1,4)  =2. 0*DELST/UI 
A (2, 1)  =UT**2/(KAP*DELT) 

A (2,2)  =0T 

A (2,3)  =UT/KAP*UI-UB 
A (2, 4) =-UT 

A (3,1) =1 . 0-OELST/DELT 
A (3,2) =-0.5*DELT/UI 
A (3,3)  =-DKU 

A (3,  4)  = (BELT  - DEI. ST ♦ DEL T*  (VT*0. 5*VB)  )/0T 
A (4,1)  =A  (3,1) -1.0 
A (4,2) =»  (3,2) 

A (4,  3)  =-DKU 

A (4,4) =DE  LST/UI*  0 . 5 *WEPF/U V FF 
COBB?  D=THETA/(XC-X) 

IF(OT.LE.O.O)  CORE  3D=0 . 0 
B (1)  = (KAP*VT)  **2*CORR3D 

B (2)  =0.0 

B (3)  = 10.0*TAUM/UI**2 

B (4) =-DWEDX-0.5*HEFP*DUEDX/UFFF 


C 

r 

C IF  UT  BECOMES  SMALL,  THEN  THE  MATRIX  BECOMES  ILL-CONDITIONED. 

C PREEZF  DUTDX,  AND  REMOVE  THE  DS  F EQN . SOLVE  THE  REDUCED  SET. 


IF  (ABS  (UT)  ,GT.0.025)GO  TO  200 
A (1,3) =A (1,4) 

A (2,  1)  =A  ( 3,  1) 

A (2,2) =A  (3,2) 

A ( 2,  3)  =A  ( 3,  4) 

A (3, 1 ) =A  (4,1) 

A (3,  2)  =A  (4, 2) 

A (3,3)  =A  (4,4) 

B (2)  =B  (3) 

B (3)  =B  (4) 

CALI.  FACTOR  ( A,  B,  IP  IVOT  , B,  3 , IFLAG) 
CALL  SnBS  T(W,B, RATES, TPIV  0‘r,  3) 
RATES  (4)  =R  ATFS  (3) 

RATES  (3)  =BUTDX 
RETURN 


2CO  CALL  FACTOR(A,H,TPI VOT ,D, 4, IFLAG) 

CALL  SUBST(B,B,RATES,IPIV0T,4) 

PFTURN 

FNB 

SUBROUTINE  DIFF2D 

CALCULATION  OF  DIFFUSERS  BITH  2-D  CORE. 

NS’’’  = 0 TS  STANDARD  DTFFUSFR,  NST=1  IS  NONSTANDARD.  SB  (J)  , J=  1 , NRC  P 1 
IS  LOBEP  BALL  VALUES  FOR  SB.  SBD ( J) , J* 1 , NHNLC  IS  UPPER  BALL  VALUES 
POP  SHU.  SHT(J)  CONTAINS  VALUES  FOR  I.OHER  BALL,  AND  J=1,NBCP1  OR 
1 , N M N LC  , BHICHF VFR  IS  LARGER. 
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PEAL*8  A (86 ,87) 

COMPLEX  C (91) 

REA  I.  X r-  (1  2 0 ) , YO ( 1 20 ) , SWT  EM P (90)  , LNV  (90) 

COMMON /BT.  IV/HS  , DFLSTS  , HU, DELSTO 

COMMON/DpR2/DELT,OB,UT,U  I , VT , V B,  UDU  I,  TAUM  , H , TH  ET  A , CELST.CFD2, 
*V  T SCO  8,  N BL 

COMMON/EFCVAL/XC  (90)  , TC  ( 90)  , A I.  (9  0)  , LH  V 

~ CMMO  N/GEOM1 /XD , W 1 ,TH, WIDTH, XI, B1, S INTH , COSTH, T WOT  HR  ,TWOTH1, 

■AS  INTH  1,  COST  HI  , AS,WH,XC1,X2,XMAX,XDE 
rOMMON/OFOM2/N,NR, NL, NO, NN1 , NLC, NRC 
COM«nN/TNITAL/DELST1, Hl,nT1,TPR1 
C OMMON/NS  TD/I C 1 , T D2 , 10 3 , NST, SWT (90) 

C OMMON/OD E1S/JSrRTS,.TENDS,NDIM,SW(90)  ,WI  (90)  , DWI  (90)  ,DDWI  (90)  , 

*08  (90) 

COMMON/ODE10/JSTRTU,JENDU,SWU (9°) ,WIU (90) 

r'OHMON/on  P2S/JT  BL  S , STBLS  (90)  , DSTAR S (90)  ,0I1DS  (90)  , BELTS  (90)  , 
SOT2DS  (9") 

COMMON /ODE2UAITBLO,STBLU  (90)  , DST  ARU  (9  0)  , 01 1 DO  (90)  ,0EITU(90)  , 
*•01200  (Q0) 

CCMMON/T.  INEAR/WDIP(°0)  ,DO2D(90)  ,DD02D  (90)  , WHD , D WE DX , UEPF, DUEDX 

COMMON /P  FSL 1/C, A ,X0,Y0 

C OMMON/PPT  NT/TPP,NOPMPR,CPEROR, ITM AX 

COMMON/SI  AD/IT, ER , CP ER R , OH  EG A 

COMMON /SPLYN/XX,WT,DWT,DDWT, ISETOP, KM  ID 

COMMON/TEMP 1/XC2, XWALLV 

COMMON/TFMP2/IEXIT,VBL(90) 

COMMON /W ALVAL/XW  (90)  , YW  (90)  ,ALW  (90) 

N LCP 1 = NLC ♦ 1 
N RCP 1 =NPC  + 1 
N MNLC  =N-NLC 


C SIMULTANEOUS  ITERATION  ON  THE  F NT I RE  CHANNEL. 

XY=0. 9 

IF  (CPE  FOR  .LE.O)  CPFROR=  .02 
TF(ITMAX. LE. 0 ) ITM AX= 1 
TT=0 

C STORE  WIDTH  OF  ENTIRE  DIFFUSER  IN  WDIF. 

DO  S J=1 , NRC PI 
S W DTF  (.1)  =WT  (J) 

C 

C CALCULATE  POTENTIAL  FLOW  IN  DIFFUSER  WITH  BARE  WALLS. 

DO  10  J=1 , N 
X C(J)  =XW  (J) 

IC  (J)  =YW  (J) 


10  A L (J)  = AI,W  (J) 

CALL  PFSL 
DO  SO  .T=NLC,NN1 

DSTARU  (.1)  =DEL  STU  *0  .00  4*SWU  ( N-  J) 

0I2DO  (.1)  = OI1*EXP(I,NV  (J)  ) 

Oil  DO  (J)  = UI2DU  (J) 

WTO(N-.I)  =UT2DU  (J) 

SO  DELTU  (J)  =0.0 

C SIMULTANEOUS  ITERATION  ON  ENTIRE  CHANNEL. 

7 W ALL  V=0 
NBL=2 

CALL  TBLSTL(O) 

CALL  CONVRT(O) 

WRITE  (6,1212) 

W RTTE  (6,  1 919) 

WRITE  (6,1111)  (K, SB  (K)  ,WI (K) ,K*1,NRCP1) 

r 

C SET  OP  THE  BOUNDARIES  OF  THE  EEC  AND  SOLVE  2-D  POTENTIAL  fLOi. 
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C**«*********nATN  LOOP  BEGINS  HERE******** 

100  IT*IT ♦ 1 

CALL  EFGEOH 
CALL  PFSL 

THE  VECTORS  X0  TO  ALONG  WHICH  THE  OI2DS  AND  OI2DO  AHE 

TO  BE  FOUND  ARE  SET  UP  BELOW. 

DO  ISO  J = 1,N 
VHAG=EXP(LNV(J)) 

ISO  X0(J)  =VHAG*UI1 
DO  250  J=1,NRC 
2S0  UI2DS  (J)  =X0(J) 

DO  260  J*NLC,NH1 
260  UI2DU  (J)  =X0  (J) 

UI2DS  (N)  =X0  (N) 

WRITE  (6,  2222) 

WRITE  (6,2221) 

COUPARE  THE  1-D  AND  2-D  VELOCITIES  AND  CP'S  ALONG  THE  Y = DELSTA B LINE. 

C PERR  = 0 . 0 
ER=0. 0 

DO  280  J=  7 , NRC 
E RR=U  I IDS  (J)  -UI2DS  (J) 

CP 1D= 1.0- (UI1DS (J)/UI1) **2 
CP2D  = 1 .0-  (UI2DS{J)/UI1)**2 
DCP=CP1D-CP2D 

CPERR= AH  AX  1 (CPERR  , ABS  (DCP)  ) 

WRTTE  (6,3111)  J ,UI 1 DS ( J)  ,UI2DS(J) , E FR, CPI D, CP2D, DCP 
280  ER=AHAX1 (ER, ABS (ERR)  ) 

DO  2 85  J=NLC, NH1 
ERR=  U1 1 DU  (J) -UI2DU (J) 

CP10  = 1 .0-  (UI1DU(J)/UI1)**2 
CP2D= 1.0- (UI2DU (J) /UI1)**2 
OCP=C  PI D-CP2  D 

C PERR=  AH  AX  1 (CPERR, ABS (DCP) ) 

WRITE  (6,3113)  J,IJI1DU(J)  ,UI2DU  (J)  , ERR,  CP  1 D,CP2D,  DCP 
285  ER=AH AX1 (ER, ABS (ERR) ) 

E RF=UI IDS  (N) -UI2DS (N) 

C P1D=  1.0-  (UT1DS  (N)/UI1)**2 
CP2D=  1,0-  (UI2DS  (N)/UI1)  **2 
DCP=CP1D-CP2D 

C PFRR=  AN  AT  1 (CPERR, ABS (DCP) ) 

WRITE  (6,3313)  N,UT1DS(N) ,0I2DS(N)  ,ERR, CP  1 D,CP2D, DCP 
ER=AHAX1 (FR, ABS (ERR)  ) 

WRITE  (6,050) ER, CPERR 


I F (CP  EPR.LE.CPEROR)  RETURN 
WRTTE  (6,1414)  IT 
IF(TT.LE.ITHAX) GO  TO  200 
WRITE  (6, 440) TT 
RETURN 

STORE  UI2DU  IN  WIU  AFTER  FLIPPING  INDICES. 

20"  DO  100  J=1 ,NHNLC 
300  WIU  (.7)  =UI2DU(N-J) 

UT  =U1 1 
T PR=I PR1 
NBL=1 
XX=0. 0 

IF  (HOD  (IT, 2)  . EO.0)GO  TO  450 
LOWER  WATL  B.L.  CALCULATION. 


U56 


U V u 


k i 


WRIT!?  (6,970) 

WRITER,  975) 

DO  120  J = 2 , NRC  PI 
120  WI  (J)  = WDIF(J) -DSTARU  (N-J) 

WI  (1)  = WDIF{1)  -DSTARU  (NM1) 

WRITF(6,  5555)  (J,SW(J)  , HI  (J)  ,WIU(J)  ,J=1,NRCP1) 

D ELST  =PFLSTS 
H=HS 

CALCULATE  PL  WITH  PRESCRIBED  PRESSURE  GRADIENT. IP  RETURNED 

VALUE  OF  IEX IT=0 , THEN  BL  HAS  REACHED  POIRT  OF  INTERMITTENT 
SEPARATION (H>=HSEP) , AND  THE  REST  HAS  TO  BE  CALCULATED  WITH 
r STRE  AMTUBF  ITER. 

DO  4?0  .1=1, NRC 
420  VEL  (.1+1)  =UI2DS  (J) 

VFL  (1 ) =UT 2 DS  (N) 

T V ALL V =0 
CALL  TPLPS 

TPfIFXIT.  NE.  0)  GO  TO  430 
WRTTF  (6,96s) 

CALL  TBLSIL(O) 

410  CALL  CONVPT(O) 

GO  TO  100 


C 

C UPPER  WALL  . USE  BL  CALCULATION  HITH  SPECIFIED 

C PRESSURE  GRADIENT ( FROM  UT2DU  OBTAINED  IN  LAST  ITERATION). 


450  T W ALL V =1 

WRITE  (6,960) 

W BITE  (6,975) 

D ELST  = DELSTU 
H =HU 

J FNDU  = NNNLC 
CALL  TPLPS 
CALL  CONVRT(I) 

GO  TO  100 

940  FOR*  AT  (*  ♦♦♦UN  ABLE  TO  CONVERGE  IN', 12,'  ITERATIONS'//) 

9S"  FORMAT  ('OLAPGEST  ABSOLUTE  ERRORS,  V ELOCITY  *• , IP  El  2 .5  (FT/SEC ) ' , 
SSX,,CPERR=',1PF12. 5//) 

960  pOPNATC  ♦♦♦♦♦♦UPPER  HALL  VALUES******'//) 

965  FORMAT  ( • CONTINUE  B.L.  CALCULATION  HITH  LINEAR  V.  P.  METHOD') 

97 0 ♦OPNATC  ♦♦♦♦♦♦LONER  HALL  VALUES  ♦♦♦♦*♦•//) 

975  FORM  AT  ( ' BOUNDARY  LATER  CALCULATION-PRESCRIBED  PRESSURE  GRADIENT') 
1111  F ORB  AT ( ' • ,15, 1P2E15.5) 

121?  FORMAT ('IDIPPUSER  WIDTH  FOR  THE  FIRST  ITERATION') 

1111  POPMAT(,-',T4,,K',T10,'SW(K)',T25,'HI  (K)  '/) 

1414  E ORN AT  ( ' 1 ♦♦♦* ♦♦ITER AT  ION  NUMBER  ',14,'  ♦*♦♦♦*#//) 

2?2?  FOPMATC  1 VELOCITY  COMPARISON'//) 

2221  FORI1AT(,0NODE#',T12,'1-D  VEL • ,T2 7, • 2-D  VEL  • ,T42  , • ( 1 D-2D  VEL)', 
«T57,'CP  1-D',T72,'CP  2-D • , T87 , • ( ID -2D  CP)'/) 

7311  FORE  AT  (•  ' , 15, 1P6E15 . 5) 

5655  FORMAT  (I5,1P1E15. 5) 

END 

SUBROUTINE  NSTAND 


r READS  IN  AND  PROCESSES  GEOHETRY  FOE  A NON-STAN  DARN  DUCT.1 

C ALSO  SUPPLIES  AN  ESTIMATE  FOR  DUCT  WIDTH  TO  BE  USED  FOR 

C SIMULTANEOUS  ITERATION  IN  THE  FIRST  LOOP. 


C OHMON/BLI V/HS ,DBLSTS, HU , DELSTU 
C OHHON/DER1/DCEL  DX, DUBDX ,DUTDX,DDIDX 

COHHON/DER2/DELT, UB,OT,UI  ,VT,VB,UDUI,TAOH,H, TH ETA , DiLST, CFD2 , 

IV ISCOS , NBL 

COMMON/EFCVAL/X  (90)  , Y (90)  , AL  (90)  , ALNV  (90) 
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no  no 


COHHON/GEOH1/XD,B1 ,TH , BI DTH, X 1 , B 1, S IBTH, COSTH, TWOTHR, TW0TH1 , 
SSIBTH1,C0STH1 ,AS,BH,XCl,X2,XHAX,XDE 
C0HH0N/GE0H2/N  , HR, BL , HO ,NH 1 ,SLC,BBC 
COHHON/NSTD/I D1 ,ID2,ID3,BST,SWT(90) 

COHHON/ODE1S/JSTRTS, JENDS,  NDIH,SW (90)  ,1/1(90)  ,001(90)  ,DDVZ( 90), 

*DS  (90) 

CONHON/ODE1 U/JSTRTU, J EBDU, SWO (90) , BID (90) 

COBH OB/IB ITAL/DELST 1,H1,0I1,IP?1 
C OHBOB/PRIBT/IPR ,BORHPR,CPEROH, ITBAX 
COHHON/S PLYB/XX, BT , DWT , DDHT, ISETOP , KHID 
COHBOB/TE  HP1/XC, I W ALLY 
COHRON/W AL V AL/XB (90) ,TB(90) ,Al«  (90) 

BRITE  (6,997) 

997  FOFHAT (' 1 ****BOB -STAB  OAR  D DOCT,  USER  IHPUTTED  BALL  CCORDIB AT  IS'// 

*•  NODE#* ,T10, »XB» ,T25, ‘TB* ,T40, » ALB'/) 

READ(5,900)N,NR,NL 

READ(  5,  91  0)  (X  B (J)  , TB  (J)  ,ALB(J)  ,J=1,H) 

WRITE  (6,9  99)  (J,XW (J)  ,TB( J)  , ALB (J) , J=1 ,B) 

IE(XB  (B)  . EQ.0.0.  AHD.TB  (H)  .EQ.O.O)GO  TO  50 
BRITE  (6,915)  XB(B)  , TB  ( B) 

STOP 

50  R BAD (5,9 10) Wl,TBOTHD, AS 
READ (5,911) B1  ,UI , VISCOS,XC 
TE(XC.EC).O.O)  XC=1.  E5 
READ(5,920)IPR,NORHPR,ITBAX,CPEROP 
READ(5,925) HS , DELSTS , HU, DELST U 
H =HS 
H 1=HS 

DELST=DELSTS 

DEIST 1=DELSTS 

IE  (HU.  BE.  0.0)  GO  TO  60 

HU=H1 

DELST  U=  DELST 1 
60  COBTIBUE 

B RITE  (6,  9 90)  B , BF , BL 

BRITE  (6,935) B1 ,TWOTHD, AS,XC 

WRITE  (6, 9 40)  B1,UI,TISCOS 

WRITE  (6,945) HS, DELSTS, HU, DELSTU 

BRITE  (6,  950)  I PR, NORNPR , CP EROR , ITHA X 

IF  THIS  IS  AB  INVISCID  CALCU  LATIOH (HOBL) , RETURB  TO  HAIH, AFTER  SETTING 

EFCV ALU ES=B ALLVALUES. 

IF  (ID 2. BE.  2)  GO  TO  70 
DO  65  J=1  ,N 
X (J)  =XW  (J ) 

T (J)  =TB  (J) 

65  AL(.T)  * ALB  (J) 

RETURN 

70  T HET  A=  DEL  ST/H 
NH1=N-1 
NRC=MR 
NLC=NRC* 1 
B RCP 1 = NRC* 1 
NHNLC=N-BLC 
H 1 = H 
UT1=UI 

DELST  1=0  EL  ST 
I PR1 =T  PR 

FTBD  ARC  LENGTHS  BETWEEN  INPUTTED  WALL  COORDINATES  USING  A 

STRAIGHT  TINE  APPROXIH ATION . 

S W (1 ) =0.0 

S B (?)  = SOP  T ( (XB  (1 ) -XV  (N))  **2»  (TB  ( 1)  -TB  (B)  ) **2) 
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no  100  J = 2 , N R c 
.1 11  =.1-1 

SW  ( .7  ♦ 1 ) =SH  (J)  ♦SORT  ( (X  W (J)-XW(JHI))  **2  ♦ (YW  (J)  -YW  (JHl)  ) **2) 

swn (i ) =n. o 

no  110  J=2, NHNLC 

,1N1  = J-1 

NMJ=N-,7 

N P=  N H J ♦ 1 

SHIJ(J)  =SB0  (JH  1 ) ♦ SORT  ( (XW (NHJ)  -XW  (NP)  ) **2  + (YW  (HNJ)  -YW  (HP)  )**2) 
SCALE=SW  (NRCP1)/SWU  (NHNLC) 

TF(NRCP1.RT.  NHNLC)GO  TO  300 

--GREATER  NO  OF  SEGS  ON  OPPER  W ALL(NRCP1  .LE. NHNLC)  . SO  HILL 
NAP  OPPFR  WAIL  COORDINATES  ONTO  THE  LOWEB  ONE. 

DO  210  J=1 ,NNNLC 
SWT  (J)  =SW0  (J)  *SC  AL  E 
WT  (1) =W1 
K = 3 

DO  259  J = 2, NHNLC 
.7  H 1=.7  - 1 
NH.J=N-,7 

TF(SWT  (,7)  . GT.  SW(  2)  ) GO  TO  230 
■-SWT(J)  IS  BETWEEN  NODES  N AND  1. 

RATTO  = SWT  (J>/SW  (2) 

X 'r  = R AT  10  * X U 1) 

Y T=R  A'FTO*  Y W (1  ) 

RO  TO  ?50 

■-SWT(J)  LIES  BETWEEN  NODE  1 AND  NRC. 

K = K + 1 

TE(SWT(J)  .GT.SW(K))GO  TO  220 
KP1  = K ♦ 1 
K H 1 =K  - 1 
KH2=K-2 

RATTO=  fSWT  (J)  -SW  (KH  1)  ) / (S  W (K)  - SW  (KH  1)  ) 

X T=XW (KN2) ♦ PATIO* (XW (KH 1 ) -XW (KH2) ) 

Y T=YW  (KH  2) +RATIO* (YW (KH1 )-YW  (KH 2) ) 

WT  (.1)  =S0RT  ( (XW  (NHJ)  -XT)  **2+  (YW(NHJ)-YT)  **2) 

WRITE (*, 955) 

WRITE  (6, 800)  (J,SWU(J)  , SWT  (J)  ,WI(J)  ,J=1, NHNLC) 

RO  TO  '"'A 

■■GREATER  NO  OF  SERS  ON  LOWER  W ALL( N RC PI . GT. NHNLC) . HAP  LOWEB 
WALL  ONTO  OPPER  ONE. 
no  Tin  ,1=1 , npc pi 

SWT  (,7)  =SW  (.1)  /SCALE 
WT  (1)  =W1 
K = 1 

no  350  J=2 . NRCP1 
GO  TO  330 
K =K+  1 

T*(SWT  (J)  .GT.SWU  (K)  ) GO  TO  320 
KH1=K- 1 

R AT TO=  (SWT  (.7)  -SWO  (KH  1)  ) / (SWU  (K)  -SWO  (KHl)  ) 

N HK  = N -K 
J H 1= J - 1 
N H = NHK* 1 

XT  = XW  (NH)  *R  AT  TO*  (XW  (NHK)  -XW  (NH)  ) 

YT=YW  (NN)  ♦RATIO*(YW(NHK)  -YW(NH)  ) 

WTfJ)  =S0RT((XH(JH1)-XT)**2*(YN(JH1)  -TT)  **2) 

DO  360  J=  1.NRCP1 
SWT  (.7)  =SW  (J) 

W RITE  (6,  960) 


n o n on 


-k- 


WRITE  (6,800)  (J,SW (J)  ,SWT  { J)  , N I ( J)  , J = 1 , N RCP 1 ) 

•>00  CONT IS 08 
J STRTS=1 
JSTRTU=1 
J ENDS=NRCP 1 
JENDU=NMHLC 

800  FORMATC  ' , 15,  1P3E15. 5) 

898  FORMAT  ( ' • ,15 , 1P3E15 . 5) 

900  FOP"  AT  (3110) 

910  FORM  AT  ( 3F  1 0 . 0) 

911  FORMA  T (4E  1 0.  0 ) 

915  FORMAT  (•  0***ORIGIN  IMPROPERLY  LOCATED, XW  (N)  = • , 1 P E 1 2.  5 , • Y W (H)  = • , 

$1  PEI  2. 6///) 

920  FORMAT  (31 10, E10.0) 

925  FORMAT(4E10.O) 

930  FORM  AT  (' 1 SEGMENT  COUNT,  TOTAL=',I2,'  LOWER  W ALL=  ',12, 

S'  UPPER  W ALL  = ',12//) 

935  FORMAT  ( *01 NLET  WIDTH=  ' , 1 PEI  2 .5, ' ( FT)  , TWOTHETA=  ,,1PE12.5, 

S'  (DEG)  , ASPECT  RATIO*  ,,1PE12.5,*  XC= • , 1 P El  2 . 5//) 

940  FORMAT('OINLET  BLOCKAGE*  • ,1PE12.5, ',  INLET  CORE  VBLOCITY=  ', 
S1PE12.5,*  (FT/SEC)  , KINEMATIC  VISCOSITY*', 1PE12.5,*  (PT2/SEC) '//) 
945  FORMAT  (• OINLET  BL  VALUES:LOWEP  H ALL-H= • , F5. 2,  • , DELSTS** , 

?>?12 . 5 , • (FT)  • /T17  , • OP  PER  W A LL-H='  ,F5.2,',  DELSTU*  •, 

$£12. 5, ' (FT) •//) 

950  FORMAT ('0B.L. PRINT  INTERVAL* ', 12,  • , PRINT  TYPE (NORBPR)  = • ,14, 

J ' , MAX  CP  ERROR*' ,F7 .5,'  , MAX  # ITER ATIO NS  = • , 1 2//) 

955  FORMAT  ('-BOUNDARY  WIDTH  FOR  THE  FIRST  ITERATION--'// 

S'  #* ,T  10,  'SWU  (J)  • ,T25, 'SWT  (J)  • ,T40 , 'HI  (J)  •/) 

960  PORMAT ('-BOUNDARY  WIDTH  FOR  THE  FIRST  ITERATION  — • // 

S*  #» ,T10, 'SW(J) • , T25 , * SWT(J) • ,T40, • WI  (J)  •/) 

RETURN 

END 

SUFROUTINF  TBLPS 

CALCULATE  TURBULENT  BOUNDARY  LAYER  PARAMETFRS  WITH  SPECIFIED 

PRESSURE  GRADIENT. 

EXTERNAL  DERPS 

REAL  KAP,  VALSM1  (3)  ,RATEM  1 ( 3) 

REAL  X P ( 3)  , YP  (3)  ,ZP  (3) 

INTEGER  KMTD, JSTART, JEND,JTBL, IRUNGE, IEXIT 
COMMON/ADAH  1/X  , V A LS  (4  ) , RA  TES  (4 ) , R ATE  (4,8) 
COMMON/BLIV/HS,DELSTS,HU,DELSTU 
COMMON/DER  1/DDEI.DX,  DUBDX  , DUTDX,  DU  IPX  1 

COMMON/DFR2/DELT,UB,UT,UI1 ,VT, VB.UDUI ,T A UM , H , TH ET A , DELST , CFD2 , 
SVTSCOS , N BL 

COMMONA  AG  AO,  TAUML  ,TAUM  EO 

COMHON/O DE1S/JSTRTS,JFNDS, NDTM.SW (9  0)  , VI  (90)  ,uVI(90)  ,DDVI(90)  , 

JDS  (90) 

COMMON/ODF 1 U/JSTRT  U, J EN  D U , SWU (90)  , VTU  (90) 

COHMO  N/ODE  2S/.1TBLS  , ST BLS  (90)  ,BSTAPS  (9  0)  ,U II  DS  (90 ) , DELTS  (90)  , 
SUI2DS  (QO) 

COMMON/ODE 2U/JT BL U, STBLU (90)  ,DSTAPU  (90)  ,UI1D0  (90)  ,DELT0(90)  , 
SUI2PU  (90) 

COM MON/PR  TNT/TPR , NORM  PR, UERR, ITM AX 
C OMM  ON/S  PI.  Y N / X X , U I ,DUIOX  , DDU  I , IS  FTU  P,  KMT  D 
COMMON/TEMP 1/XCMX, IWALLV 
C0MMON/TE"P2/TEXTT,VEL (90) 

*•0017  Af  FNCE  (RAT  EM  1 ( 1)  , DDELDX)  , (VALSM  1 (1)  , DELT) 

TMPOSFD  PRESSURE  GRAD  IS  OBTAINED  VIA  CALL  TO  SPLINF. 

THE  FOB  T V AT.  ENCTN  0 IMPLICITLY  SETS  DELT, OB, UT  EQUAL  TO  VALS. 

SET  UP  THE  EXTERNALLY  IMPOSED  PRESSURE  FIELD 

IP  (IW  AT.LV  . EO.  1)  GO  TO  20 


TP  (XX.  IT.  0.0)  XX=  S W (1) 

X=XX 

▼1AX  = SW  (J  PHDS) 

I SPTUP=0 
KSID=JSTRTS»1 

CULT.  SPLINE  (SW,VEL,DVI,DDVI,DS,  JST  RTS, J EN DS, N DIB,  4) 

RO  TO  TO 

T P (X  X • T.T.  0. 0)  XX=SWU(1) 

X=XX 

X NAX=SWU  ( J END  II) 

t setups 

KNID=JSTRTU+1 

CALL  SPLINE  (SHU,  VTO  ,DV I, DDV I , DS, JSTRTD,  J ENDU, N DIH , 4) 


INITIALIZE  COUNTERS,  CONPUTE  START  VALUES  OP  UT, UB, CELT, 

TO  ,T T BT.  = 0 
NLOOP=0 
U IREP  =UT 
CALL  START 
CALL  TAUNAX(TAUN) 

»n=x 

HO=H 

T AUNL=TAUHEQ 
’MUNsTAnNEQ 
U RITE  (6,  900) 

IRUNGE=1 
DX=DELST 
I EXIT=° 

DO  SO  J = 1,  T 
VALS ( J) =VALSH1 (J) 

R ATEN  1 (J)  =0.0 
SO  OATES  (.1)  =0.0 

H SEP=  1 . ♦ 1 . / (1 . -DELST/DELT) 


C npGIN  HATN  LOOP. 

c PRINT  INITIAL  VALUES. 

GO  TO  1nS 
100  NL00P=NL00P»1 

H SEP= 1 . ♦ 1 . / (1 . -DELST/DELT) 

c IP  IWALLV=0 (LOWER  BALL)  AND  H>=0.9*HSEP  OR  H<HO,S»ITCH  OVER  TO 

C LINEAR  V.P.  ITERATION.  RETURN  TO  CALLING  ROUTINE  APTER 

C SEATING  B.L.  VALUES  TO  THAT  AT  THE  LAST  PRINTOUT. 

IP  (IWALLV.EO. 1)G0  TO  106 

CP=1 . 0-  (UI/UIREF)  **2 

TP  (CP.LE.0.2) GO  TO  107 

IP  (H. LT.O ,9*HSEP. AND. H.GT. HO) GO  TO  107 

WRITE  (6,920) X 

920  *ORfl AT  (•  ***♦ H>. 9*  HS  EP  OR  H<HO  AT  X=*,1PE12.5) 

H =H0 

UB-UBO 

(1T=UT0 

Dumxi=nuiDXo 

TAUH=TAUHO 

TPfTWALLV.EQ. 1)  GO  TO  104 
XX=S™BLS ( JTBL) 

DELST=DST ARS (JTBL) 

D ELT  = DELTS  (JTBL) 

U I1  = U T1DS  (JTBL) 

JTBLS=JTBL 

RETURN 


nodf*  or>oo  no 


109  XX-STBLO  (JTBL) 

DBLST»DSTAR0  (JTBL) 

DBLT=  DBLTO  (JTBL) 

0I1-0I1D0  (JTBL) 

JTBL0*JTBL 

RBTORR 

IF  H>=HSEP,  AND  OPPBR  MALL,  THBR  SET  DFLST*CONST  TO  EXIT. 

106  IF(H.tT. 0. 9*HSEP) GO  TO  107 
WRITE  (6,030)  X 

930  FORBAT ('  IRTERBITTERT  SEPAR ATIOR  AT  X=«,1PE12.S, 

S'FT,  SET  DELST=CORST  TO  EXIT'/) 

JTBLO*JTBL 

RETORR 


STORE  CORRERT  VALDES. 

107  DO  110  J=  1 , 3 

V ALSB  1 (J)  =VALS  (J) 

110  RATEB1  (J) =RATES  (J) 

STORE  THE  LAG  PARAHETERS. 

XO=X 

T A0BL=T AOR 

PRINT  CORRENT  VALUES. 

XX=X 

IF(BOD  (NLOOP,  IPR)  ,NE.O)GO  TO  150 
105  JTBL= JTBL* 1 

HS  EP= 1 . ♦ 1 . / ( 1 . -DELST/DELT) 

DQDX=10,0*TAOB/OI**2 
CP=1 . 0- (DT/DI REF) **2 

WRITE  (6,910)  XX, DEL  ST,  H,H  SEP,  CP,  DELT,  0 B,  DT,  01,  CFC2,DQDX 
IF (IW ALLY . EQ. 1) GO  TO  120 
STBLS ( JTB  L)  =XX 
DSTARS  (JTBL)  = DELST 
DELTS  (JTBL)  = DELT 
0I1DS  (JTBL)  =01 
GO  TO  130 
120  STBLD(JTBL)=XX 

DSTARO  (JTBL)  = DELST 
DELTO  (JTBL)  = DELT 
0I1D0  (JTBL) =01 
130  HO=R 
tlBO=nB 
nTO=OT 

DOTDXO=DOIDX 
T AOHO=TAUB 


IF (IEX IT . EQ. 1) GO  TO  260 
150  CALL  A DABS  (DX  , 3 , DERPS , IRORGE) 
IP(X.  LT.  XB AX) GO  TO  100 


EXIT  VALOR  COBPOTATIORS. 

OBTAIN  VALOES  AT  X=XBAX  BY  EXTRAPOLATION. 
I EXIT* 1 
XY=XH  AX 

IFdWALLV.EQ.  1)GO  TO  165 
DO  160  J= 1 , 3 
J J=JT  BL- 3 *J 
XP(J)  = STB  t S (J  J) 

T P (J)  = DSTARS  (JJ) 

160  ZP(J)  *01 1 DS  (J  J) 
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r,o  Tn  170 
165  00  16  8 .1*1,3 
J J=.1TBL-3*J 
XP(J)  =ST BLO  (JJ) 

YP  (.1)  =DSTARO  (JJ) 

168  ZP(J)  =OT  1 DO  (J  J) 

170  DELST  = YINT  (XP, YP,XHAX) 

U I=OT2PO  (JENDO) 

JTBLS=J’’'BL*1 
.IT  BLO= JTBL  ♦ 1 
GO  TO  106 
260  CONTINOF 

000  PORHATflHO, • X DSTAR  H HSEP  CP  DB 

OB  OT  01  CP/2  DQDX/OI'/) 

0 10  POP1AT{2F10.5,4X, 2F7.3,2X,P6. 3, 4F 1 2 . 5 ,P 1 2 .6 , Pi 2 . 5) 

R FTflR  N 
FID 

SOBPOOTINE  BLVALO 

C GIVEN  OT.OT, OB, DELT,  COMPOTE  B . L.  THICKNESSES  DSTAR  ,D2STAK, 

C AND  SHAPE  PACTORS  H AND  HBAR.  ALSO  CP/2*CPD2. 

COHNO  N/3L0LD/ET0LD ,01 INO,OZOLD 
COBBON/DER1/DDELDX,DOBDX,DOTDX,DOIDX1 

C0HH0N/DPR2/DFLT, OB, OT ,011 ,VT, VB,ODOI,TAOH, H,  THETA, DSTAR, CPD2, 
SVTSCOS.NBL 

COHBON/L AG/XO ,TAOHO,TAOHEQ 
COHHON/SPLYN/X,OI ,DOI,DDOI,ISETOP,KHID 
REAL  PI,  P ID2 

DATA  PI, PID2/3. 141503, 1.57079/ 

VT=0T/(0. 41*01) 

VB=OB/OT 

DSDDEL=VT»0.5*VB 
D ST A R=DSDDEL*DELT 

THDDEL=DSDDEL-2. 0*TT**2-0. 375*VB** 2-1 . 50949* VT* VB 
THETA=THDDEL*DELT 
H=DS'rA  R/THETA 
CPD2=  (OT/OI)  **2 
IF  (OT.LT.O.O) CED2=-CPD2 
RETOR N 
C 

c ************* ******** 

c 

ENTRY  START 


C GIVEN  H, DSTAR, 01,  FIND  OT,  OB.  ZI*DSTAR/DELT. 

C OBTAINED  BY  PITTING  COLES'  WALL/WAKE  PROFILE  TO  THE  INPOT  . 

C INITIALIZE  LAG  E0OATION  PARAHETERS,  AND  TA01. 

XO=X 

IE(TAOH.LE.0.0)TAOH  = 10.0 
TAOHO=TAOH 

R POST*  01*  DST AR/V IS COS 

C USE  FLAT  PLATE  VALUES  FOR  FIRST  GOESS. 


VT=0 . 016B5/REDST**0. 166667 

7,1=0.125 

VB= (ZT-VT) *2. 0 

HDH11 *H/  (H-1 .0) 

NOOT=0 
00  NL2=0 
NLOOP=0 

A DST=OST  AR*0. 4 1*OI/VISCOS 
100  VTO*VT 
Z TO=Z  I 

NL0OP=NL00P*1 
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TF(10.GF. NLOOP) GO  TO  150 

WRITE  (6,900)H,DSTAR,0I,0T,0B, VT 

STOP 

ISO  FVT=VT* (2 . 05*  ALOG (ADST*ABS  (VT)/ZI) ) *2.0*(ZI-VT) -1.0 
FPVT=1.0S*ALOG  (ADST*  ABS  ( VT)  /Z I) 

VT=VT-PVT/PPVT 

IP(0. 0001 .LT. ABS (1.0- VTO/VT) ) GO  TO  100 

C BEGIN  LOOP  FOB  VB  ITERATION. 

170  V PO= V B 

NL2=NL2*1 

IF(10.GE. NL2) GO  TO  1B0 

WRITE  (6,915)  H ,DSTAR  ,tJI,OT,OB,ZT 

STOP 

180  PVB=HDHM1 * (2. 0*VT* VT* . 175* VB*VB* 1 . 5 98949* VT* VB) -VT-0.5*VB 
F PVB=HDHH 1 * (0 .75*VB*1 . 59  89  49*VT) -0.5 
VB= VB-FV8 /FPVB 

IF(0.  0001.LT. ABS (I.O-VBO/VB)  ) GO  TO  170 
Z I=HDHH1  ♦ (2 .0  *VT*VT*0 . 17 5*VB *VB* 1 . 5 98 94 9*VT* VB) 

IF(0.  0001  .GT. ABS (1 .0-ZT/ZIO)  ) GO  TO  200 
NOHT=NO(TT*  1 

IF(NOUT.LE. 10) GO  TO  99 

WRITE  (6,910)H,DSTAR,UI,OT,OB,ZI 

STOP 

200  UT=VT*.41 *01 
OB=VB*OI 
D FLT=  DST  AR/ZT 
CFD2=  (OT/OI)  **2 
TPfUT.LE.  0.0) CFD2=-CPD2 
WRITE (6,920)OT,OB,DELT,DSTAR, H 

INITIALIZE  STARTING  GOESSES  FOR  TAOMAX. 


t ^ 


ETOLO=0.  25 

900  FORMAT  (’  VT  FAILED  TO  CONVERGE  IN  10  ITERS  , H, DSTAR , 01 , OT ,OB , VT 

$ = • , 6E1 2.5) 

910  FORMAT  (•  ZT  FAILED  TO  CONVERGE  IN  10  ITERS  , H, DSTAR , U I ,OT ,UB , Z 

«T  = • , 6E 1 2 . 5//) 

915  FORM  A T ( • VB  FAILED  TO  CONVERGE  IN  10  ITERS, H, DSTAR , 0 1, OT, UB , VB=* , 

$ 6 FI 2. 5//) 

9?P  FORM  AT  (•  START  VALOES,OT=  ',F10. 5,*  0B=  •,*10.5,*  DFLTA=  •, 

SP10.5,'  DELST=  • ,F10.5,*  H= 

9FT0FN 

********** 

FNTRI  TAOMAX(TAOLAG) 

GIVEN  ODOT  AT  WHICH  TAO  IS  MAX, FIND  THE  CORRESPONDING  ETA=Y/DELT 

, MA X DODY,  AND  THE  MAX  SHEAR  T AHM/RHO. 

UDni  IS  SET= .76  FOR  ATTACHED  FLOWS  AND  =0.6  FOP  DETACHED  FLOWS. 

0001=9.76 

TF  (OT.LF.0.O)  T!Dfir=0.60 

n DM!  = 1 . 0- ODHI 

FT=ETOLO 

TF(OOT.LT.O.O) GO  TO  290 
*T=0. 25 
GO  TO  120 
290  NL=0 
100  =TO=ET 

TF(ET.GT.O.O)  GO  TO  110 
WRTTF  (6,°60)  FT 
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ET=0. 25 
RO  TO  120 
110  NL=NL>1 

T E ( NL . GT  , 10)  GO  TO  590 

E P','= VT*A  LOG  (ET)  - VB*  (COS  ( PI  D2  *ET)  ) ♦♦2+UDH1 
PPET  = VT/FT*PT  D2* VB*  SIN (PI*ET) 

ET=ET-EFT/FPET 

TP  (0. 0001  . LT. ABS (1 .P-ET/ETO)  ) GO  TO  100 

TP (ET . tT . 0 .25 ) ET=0 .25 

ETOT,D=FT 

120  DUDY=  (UI/DELT)  * ( VT/ET* VB  *PID2*S  I II  (PI*  ET)  ) 

USING  KUHN-NIELSEN'S  EDDY  VISCOSITY,  BHICH  INCLUDES  EFFECTS 

OF  TNTFRfHTT PNCY  AND  PRESSURE  GRADIENT  PARAHETER, 

BETA  IS  THE  CLAUSES  PRESSURE  GPADIENT  PARAHETER.  EPS  IS  EDDY 
C VTSCOSTT Y , GANNA  IS  INTERHITTENCY . 

C SET  EPS  = 0 .011  THF  *RBF  SHEAR  LI  HI T FOR  FLOWS  THAT 

C ARE  NEAR  AND  BEYOND  DETACHNENT. 

E2=0. 0 

TP  (UT.LF.O.S) GO  TO  140 

F2=0. 001 R 

RETA=-DSTAR*UI* DU  1/(15. 0*UT*UT) 

T F (DUI.GE. 0. O.np. ABS(BETA)  .GE.174.0)GO  TO  340 
E 2=o. 001 R*EXP (-BETA) 

140  E PS=0 . 0 1 1 *E2 

GANNA =1 .0/  (1 . 0+0 ,C*ET**6) 

T AnNEQ=P  PS*GA NHA*DUDY*UI*DSTAR 

C INCLUDE  SHPAR  STRESS  HISTORY  BY  LAGGING  THE  EQUIL  STRESS 

H LAN  = 0 . 02  5 

TF(Ut.i,p.  0.0)  H LA  H = 0 . 70 

DTNDX  =HI.  A N*  (T  AON  Ep-TAH  NO) /DELT 

T AUL  AG=T AUNO*  DTN  DX* (X- TO) 

R ETUR  N 

SO 0 WRITE (6, 950) NL,ET,UT,UI,UB 
STOP 

050  FOREST^HC,1  E^A  FAILED  TO  CONVERGE  IN' , 14, • ITERATIONS, ET, UT, UI,UB 
t=  ',4P15.5) 

oso  E ORN AT  (•  ET  SET=0.25,  OLD  VALUE  WAS  ■ ,F10.5) 

END 

SUBROUTINE  PSTEST 

C DRIVER  ROUTINE  TO  TEST  ■’’n  RBU  LE  NT  BOUNDARY  LAYER  CALCULATION 

C WT^H  SPECIFIED  PRESSURE  GRADIENT. 

r SUBROUTINES  RFOUIRED:  A DANS , BL VAL U ,DERPS , FACTOR , PSTEST, BKS4 , 

C SPLINE, SUBST ,TPIDIAG. 

C ♦♦♦♦ROUTINE  TO  TEST  NEW  BOUNDARY  LAYER  PREDICTION  HETHOD 


C 

C 

C 


USTNG  TAUNAX-ENTRAINHENT  CORRELATION. 


CONNON/BLIV/HS,DELSTS,HU,DELSTU 
COHNON/DFR  1/DCPLDX,DHBDX,DUTDX,DUIDX1 

CONNON/DEP?/DELT,UB,UT,UI1 ,VT,VB,UDUI,TAUH, H,TH ETA, DEL ST, CFD 2, 
$V  TSCO  S, N BL 

COHNON/ODE  IS/ J ST RTS ,JENDS,NDIB,SW(90)  ,VI(90)  ,DVI(90)  ,DDVI(90) , 
*DS  (40) 

COHNON/ODE 1U/JSTRTU,J END U,SWU (90) , VIU (90) 

CONHON/ODE 2H/.1TBLU  ,STBLU  (90)  , DST  ARU  (9  0)  ,U  II  DU  (90)  , DELTU  (90)  , 
SUT2D0  (40) 

COHNON/PRTNT/TPR,NOPHPR,CPFROR,ITHAX 
CONNON/SPLYN/XX, UT , DUI , DDUT , IS  ETUP , KH ID 
COHHON/TEHP1/XC,IWALLV 
NDIN*90 

— -ENT»p  THE  IMPOSED  VELOCITY  DISTRIBUTION 
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X=0.0 

READ  IN  THE  B.L.DATA.  END  I.AST  CASE  WITH  2 BLANK  CARDS. 

DELST  AND  THETA  ARB  IN  FT. 

TAUM  IS  THE  STARTING  VALUE  OF  THE  MAX  SHEAR  STRESS. 

READ(6,902)NPTS 
IE  (NPTS.LE.O)  GO  TO  800 
READ (5,90lt) XX, DELST, H.VISCOS 
READ  ( 5,  905)  IPR,XC 
IE (XC.EQ. 0. 0) XC=1 .E5 
IE (IPR.LE.O) IPR=  2 

WRITE  (6,006) XX, DELST, H, VIS  COS, XC,IPR 
R FAD  (S  ,0  1 0)  (SWU(I)  ,VIU(I)  ,I=1,NPTS) 

WRITE  (6,916) 

W RTTF  (6,  920)  (SWO(I)  , VIO(I)  ,I  = 1,KPTS1 
THETA *DELST/H 
J STRT  0=1 
.1 ENDO  =NPTS 
DELSTtT=DE  LST 
HU=H 

I WALL V=1 

0I2D0  (JEN DO)  = Vin(JENDU) 

CALL  TBLPS 
W RITE (6 , 9 10) 
non  CONTINUE 
90?  FORMAT  (1 1 0) 

9"U  POPMA  T (UE 1 o, 5) 

906  FORMAT  (IIO.EIO.S) 

906  F0PMAT(1Hr',,X, DELST, H,  VISCOS  , XC,  IPR=  » , 1 PSE 1 2.  9,  I 5///) 

910  FORMAT  (8  F 10.0) 

916  FORMAT  (1 H 0 , • X UI'/) 

920  PORN  AT  { 1H0 , 2E20.  5) 

910  PORMAT  (1HO,«***TN  MAIN  ROUTINE***') 

RETURN 
END 

SUBROUTINE  TNVCID 
CALL  PESL 

CALL  OUTINT  (TPOINT) 

RETURN 
END 

SUBROUTINE  OUTINT (IPOINT) 

C COMPUTE  FUNCTION  VALUES  AND  DERIVATIVES  AT  INTERIOR  POINTS. 

R EAL*  8 A (86,87) 

PEAL  X0(120)  , 10(120)  ,LNV( 90), TX 0(1 20)  ,TYP  (120) 

COMPLEX  FPZ  (9  0)  , ETA  (120)  ,C  (9  1)  ,FPZ0  (120)  ,FPPZ0  (12  0)  ,20  (120) 

COMPLEX  ZO,ZEPO,ITWOPT,CMPLX,CEXP,CONJG,  EPZO,  PGP  AD,  PGR  ADS 
COMPLEX* 16  GAMMA  1 (120) .GAMMA 2 (120) 

COMPLEX* 16  ET AJ , FT  A JP , ET A JM, FPP, LE RP, CDLOG , P DZO , FDDZO 
COMMON/EPCVAL/XC (90) , TC(90) , AL(90) , LNV 

COMMON/GEOM  1/XD.XI.  ,TH,WIDTH,X1,B1,  SI  NTH  ,COSTH  ,TWOTHR  ,TWOTH1  , 

SSTNTH1 .COST HI ,AS,NH,XCE,X?,XMAX,XDE 
COMMON/GF.OM2/N  ,NR,NL,NU,NM1  , NLC  , NPC 
COMMON /TNTTAL/D EL ST1, HI, VIN.IPR1 
C OMMON/PFSL1/C, A , Xn , TO 
COHMON/UPTNT/TPR,NORMPR,UERP , ITMAX 
COMMON/WALVAL/XW  (90) ,TW(90) ,ALW(90) 

Z ERO=  (0.0, 0.0) 

I TWOPI = (U  .0,6  ,281186) 

PEAD  (6,901) LINES 

C LIN"S  = *LTNES  ALONG  WHTCH  ANALYTIC  FUNCTION  AND  DERIVS  ARE  CALCULATED. 

IE  (LT  NFS . EQ.O)  PRTUPN 
IP (LI NffS. LT. 0) GO  TO  716 
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I POT  N T=° 

no  714  L*1, LINES 

ENTER  1 CURD  FOP  EACH  LIME  ALONG  WHICH  THE  INTERIOR  POINTS 

APP  TO  BK  CALCULATED.  (X  1 , Y 1 ) , (X2  , Y2)  , ARE  START  AN!  END 
POINTS  OF  LINE,  AND  NSFGS  IS  NO.  OF  SEGMENTS  ALONG  LINE. 
READ (5, 951) XI , Y1 , X2 , Y2 , NSEGS 
IPfNSEGS.  GT.O)  GO  TO  712 
I POI NT=I POINT* 1 
X 0 (IPOINT) =X1 
TO  (IPOTNT) =Y1 
GO  TO  714 

71?  NPOTNT=NSEGS+ 1 

DX=(X2-X1) /NSEGS 
DY=(Y2-Y1)/ NSEGS 
DO  71  3 >1=  1 , NPOINT 
.1 M 1 = J - 1 

T POTNT=I POINT* 1 
X 9 (IPOINT) =X1 »DX *JN1 
717  YO  (IPOINT) =Y1*DY*JN1 

714  CONTINUE 

715  N PI  = N+1 
VSCALE  = VIN 
DO  717  >1=1  ,N 

717  PPZ(.I)  = CEXP  (CNPLX  (LNV(J)  ,-AL(J))  ) 

WRITE  (6, 950)  IPOINT 
DO  760  K=  1, IPOINT 
ZO  = CNPLX  (XC  (K)  ,Y0  (K)  ) /XL 
7 0 ( K ) = ZO 

C (NPl  ) =C  (1) 

DO  720  J=  1 , HP  1 
7? " ETA  (J)  = C (J)  -ZO 
DO  710  .1  = 1, N 
ETAJ  = ETA  (J) 

'TA.IP  = ETA  (J  ♦ 1) 

ERP  = FT  A J P/ETAJ 
LERP  = CDLOG(EPP) 

GANNAI(J)  = LE  RP/ (ET AJP-ET AJ) 

G Add  A3  (.1)  = 1 . 0/  ( ET  AJ*ETA  J P) 

77"  CONTINUE 

FDZO  = ZERO 
FDD70  = ZERO 
nO  740  .1  = 1 , N 
JP1  = .1*1 
JMl  = .1-1 

IF(JHI.EQ.O)  JN1  = N 
FTA.l  = ETA  (J) 

ETAJ  P = ET  A (.7  PI) 

ETA.1N  = ETA  (J INI ) 

FDZO  = FDZO*FP7 (J) • (ETAJP*GANNA1  (J) -ETAJH*GANNA 1 (Jfll) ) 

74  n FDDZO  = FDDZOi-FPZ(J)  * (GA NNA 1 ( JN1 ) -G ANN A1  (J) ♦ ET AJP* GANNA  2 (J) - 
f>  ETA  JH*G  AHNA2  (JN1)  ) 

PPZO(K)  = FDZO/ITWOPI 
760  FPPZO(K)  = FDDZO/ITWOPI 
850  IF  (NORNPR)  860,960,870 

NORMALIZED  NEUMANN  PRINTOUT. 

860  WRITE  (6,95?) 

WRITE  (6,  960) 

DO  965  R*  1 ,1  PCI NT 
FPZO  = FPZO(R) 

V NAG  * CABS (FPZO) 

!JT  = REAL  (FPZO) 


no  n n on  non 
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VI  = -AI HAG  (FPZO) 

ALPHA  = ATAN2  (VI, UI) 

STORE  VELOCITY  COHPONENTS  LOCALLY  PARALLEL  TO  THE  HALLS 

IP  TXO  AND  TYO  APPAYS. 

AHA=ALW(K) -ALPHA 
TXO(K)  =VSAG*C05  (AHA) 

TYO  (K)  =VHAG*SIN(AHA) 

865  HPTTE  (6, 955)  K ,Z0  (K)  . 01,  VI , VHAG , ALPH  A ,K 
IE  (HORHPR)  880,870,870 

DTHEKSTONAL  NEDHANN  PRINTOUT. 

870  WRITE  (6,954) 

H PIT  E (6 , 9 6 1) 

TO  875  K==  1 , 1 POIN T 
FPZO  = FPZO  (K)  * VSCAL  E 
VHAG  = CABS (F  FZO) 

HI  * REAL  (FPZO) 

VI  = - AI H AG  (FPZO) 

ALPHA  - AT AN2  (VI  01) 

878  WRITE  (6,956)  K , XO  (K)  , YO  (1C)  , UI  , VI  , Vfl  AG  , ALPH A ,K 

PRESSURE  AND  PRESSURE  GRADIENT  CALCULATIONS. 

880  IF  (NORHPR)  882,882,885 

NORM  ALT7EC  PRESSURE  DATA. 

882  WRITE  (6,952) 

W RITE  (6,  970) 

DO  884  K=  1 ,1  PCI  NT 
FPZO  = FPZO (K ) 

VHAG  = CABS (FPZO) 

V HAGS  Q = V(1  AG*VH  AG 

CP  = (P-PTN)/0IN  = 1 -(V/VIN) **2 

CP  = t .0- VHAGSQ 
PGRAD  = -FPZO*CONJG (FPPZO (K)  ) 

PGRAOS  = PGP  A D*PPZO/VH AG 
CURVE  = -AIHAG(PGRADS) /VHAGSQ 

884  W RITE (6,  965)  K, Z 0 (K ) , PGR  AD ,PGR A D S, CUR VE ,C P, K 
IF  (NORHPR)  1000,885,885 

DIHENSIONAL  PRESSURE  DATA. 

885  WRITE  (6,884) 

WRITE  (6,975) 

VINSO  = VIN*VIN 

V SCX L = VSCAL E/XL 
DO  880  K= 1 ,1  POINT 
FPZO  = FPZO  (K) *VSCALE 
VHAG  = CABS  (FPZO) 

VHAGSQ  = VHAG  *VH AG 
PDIFF  = ( V IN SQ-VH  AGSQ ) /2 
PGRAD  = -FPZO*CONJG (FPPZO (K)  ) *VSCXL 
PGRADS  = PGRA  D*F  PZO/VH  AG 
CURVE  = -AIHAG(PGRADS) /VHAGSQ 

890  WRITE  (6,966)  K , XO  (K)  , YO ( K)  , PGRAD, PG RA DS, CUR V E, P DI FP , A 

RETURN  THE  VELOCITY  COH PON ENTS  LOCALLY  PARALLEL  TO  THE  WALL 

IN  XO  AND  TO. 

1008  DO  900  J=1  ,N 

XO(.I)  =TX  0 (J)  * VSCAL  E 

900  YO  <>T)  =TY0  (J)  *VSCALE 
RETURN 

901  POPHAT  (T  1 0) 

980  FORNATC 1HP/1HP,24X,' VALUE  OF  ANALYTIC  FUNCTION  AND  ITS  DERIVATIVES 
F.  AT ' , T4,  • BOUNDARY  AND/OR  INTERIOP  POINTS.') 

951  POPHAT  (4*10.0  , TIC) 

98?  FORHAT(1h',,61X,'  NORHALIZED  VALUES') 

954  PORH  AT  ( I HO , 67X , ' DIHFNSIONAT,  VALUES') 
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966  FORMAT  (1H  , 27  X , 1 4 , 6F  1 4 . 6 , 1 3) 

966  pn«?KAT(1H  ,27X,I4,1P6E14.5,I3) 

96  9 pops  AT  (1  HO, 3"  X,'  #•  ,PX,  *X0»  , 12X,  • YO*  , 1 2X  , 'U* , 1 3X  , • V • , 12X, 


F,<  V EL-MAG*  , 8X,  * ALPHA*  , 3X,'  **) 

461  FnRIATnHO^OX.'t'^X,  • XO  * , 1 2X , • YO  » , 1 2X  , * U • , 1 3X  , • V * , 12X  , 

R* VEL-MAG'  ,7X,  • ALPHA *, 5X ,*#• ) 

466  ° ORN  AT  ( 1H  , 13X,14,8F14.6,I3) 

966  FORMAT  (1H  , 1 3 X ,1 4 , 1 P8  E 14 . 5 , T3 ) 

4 7n  fORH AT  (1  HP, 16  X,*  #* ,9X,'X0* ,12X, • YO* ,7X, • (DP/DX)  /RHO*  ,3X,* (DP/DT)/R 
&H0’,3X,*  (DP/OS)/RHO*  ,3X,’  (DP /ON) /RHO* ,5X , ‘CURVATURE*  ,8X,*CP*  ,6X, 

S*  #•) 

976  P0R(1  AT(1H0, 16X,' #• , 7X, *X0» , 12X, ' YO* , 8X,'  (DP/DX)  /RH0*,3X, • (DE/DY)  / 
6PH0* , 3X, • (DP/DS)  /RHO* , 3X, * (DP/DN ) /RHO • , 5X , • CURV ATURE' ,3X, 

6*  (P-PIN)/RHO*  , 2X , ' #•) 

PND 

SUBROUTINE  CONVRT(IWALL) 

INTERPOLATE  AND  CHANGE  INDICES  FFOfl  JTBLS  OR  JTBLD  TO  N. 

DONE  PY  CALLING  SUBROUTINE  CHANGE. 

CONHON/CON/SVAL (90) ,YVAL(90) 

C ONHON/GEOH2/N ,NR,NL,NU,NN 1,NLC, NEC 

C0HH0N/0DE1 S/ JSTRTS  , JENDS  , NDIPt , S H ( 90)  , W I (90  ) , DU  I (90)  , DDWI  (90) , 

IDS  (90) 

C ON  HON/ODE  1U/JSTRTU,  JENDU,SWU  (90)  ,HIU  (90) 

CONMON/ODE2U/JTELU,  STBLU  (90)  , DST  ARU  (90)  ,T1IlDU  (90)  , DELTU  (90)  , 

* U I2DU  (90) 

C OMWON/ODE2S/JTBL  S , STB  LS (90 ) , DST A PS (90)  , U1 1 DS (90 ) , DELTS (90)  , 

!U  T 2D  S ( 90) 

COMNON/SPLY  N/XTNT ,FINT,EPINT, PPP INT , ISETUP , KH ID 
rOflNON/TPHPl /XCfIX  .IWALLV 
I ” ( IH  ALL. EQ. 1) GO  TO  600 
IF(.1TPLS.GT.90)gO  TO  1600 

T1=DSTARS  (1) 

•”2=DE!.TS  (1) 

T 3 = UT  IDS  ( 1) 

DO  509  K= 1 ,3 
I w (K-  2)  100,  20  C,  300 
100  DO  160  J=1, JTBLS 
S VAT.  (J)  =STPLS  (J) 

169  v VAL  (J)  =DSTAPS  (J) 

CALL  CHANGE 
DO  160  J=1,NRC 
160  D STARS  (J)  = YVAL  (J) 

GO  TO  609 

290  no  260  .7=1, JTBLS 
260  YVAL  f J)  = DELTS  (J) 

CALL  CHANGE 
no  260  J = 1 , N R C 
260  DELTS  (J)  =YVAI,  (J) 

GO  TO  500 

399  DO  960  J=1 ,JTPLS 
350  YVAL  (J)  =UI10S  (J) 

CALL  CHANGE 
no  360  J = 1,NRC 
36**  mms  (J)  = Y VAL  (J) 

609  o ONT  T NOP 

r SMARTING  VALUES  AT  NODE  N. 

9S*IBS  (N)  «T1 
Oft  TS  (N1  *T2 
• » ***S  (•)  »T3 
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RRTORN 

600  IP(JTBL0. GT.90)  GO  TO  1600 
DO  1000  K=1,7 
I? (K -2)  700,800,900 
700  DO  750  J* 1 , JTBLO 
SYAL(J)  =ST  BLO  (J) 

750  YTAL(J)=DSTAPO  (J) 

CALL  CHANGE 
DO  760  J*  NLC, NN1 
760  DSTARO(J)  =YVAL  (J) 

GO  TO  1000 
800  DO  850  J= 1 , JTBLO 
85  0 Y VAL  (J)  =DELTO  (J) 

CALL  CHANGE 
DO  860  J=NLC,NB1 
860  DELTO  (J)  =YVAL  (J) 

GO  TO  1000 
900  DO  950  J=  1,  JT  RLfJ 
950  YVAL  (J)=0I1D0  (J) 

CALI,  CHANGE 
DO  960  J = NLC, NH 1 
960  OIlDO  (J)  =Y  V AL  (J) 

100O  CONTINOE 
RETORN 

JTBLO  OR  JTBLS  IS  GT  90.  PRINT  ERROR  NESSAGE 

AND  STOP.  CORRECT  BY  INCREASING  IPR. 

1500  WRITE  (6,920) JTBLS 
STOP 

1600  WRITE (6,910) JTBLO 
STOP 

9?9  PORBAT  (• -JTRLS=' , 12, ', WHICH  IS  .GT.90,  INCREASE  IPR  AND  RERON'//) 
9ir>  PORHATC -JTELO='  ,12,'  , WHICH  IS  .GT.90,  INCREASE  IPR  AND  RERON'//) 


END 

SOBROOTINE  TBT.ST  (IWALL) 

— CALCOL ATE  DSTARS,OI1DS  FOR  A TBL  AND  1-D  CORE  IN  SIHOLTANEOOS 
ITERATION. 

TWA  LL=0  IS  LOWER  WALL.  IWALL=1  IS  OPPER. 

EXTERNAL  DEPSI 

P p AL  KAP,VALSB1  (4)  ,RATEH1  (4) 

REAL  XP  (7)  , YP  (1)  , ZP (1)  ,OP(?) 

INTEGER  KBID,  JSTART,  J END,  JTBL  , IRONGE,  IEXIT 
COBBON/ADABl/X  ,VALS  (4)  , RATES  (4)  , RATE  (4 ,8) 

CORBON/PERl/DEELDX ,DOBDX ,DOTDX,DOIDX 

COBNON/DE  P 2/DELT,  0B,0T,0I, VT, »B , ODOI , TA OB ,H ,THET A, DELST ,CFD2 , 
*YISCOS,NBL 

COHHON/GEOH 1/XD, W 1 ,T H, WIDTH, X 1 , B 1 , SINTH ,COSTH ,TWOTHR,TWOTH1, 
♦SINTH1 ,COSTH1 , AS ,WH,XC,X 2, XBAX ,X DE 
COHNON/GPOB?/N,NR , NL, NO, NHl , NLC, NPC 
COBNON/TN IT AL /DELST 1, HI, Oil, IPR 1 
CONNON/LAG/XO ,TAOBO ,TA0BEQ 

C ORB ON/ODE  IS/ JST RTS ,JENDS,NDIN,SW(90)  ,WI (90) ,DWI  (90)  ,DDWI(90)  , 
YDS  (90) 

COBB ON/O DEI O/JSTRTO ,JENDO,SWO (90) ,WI0  (90) 

C OHBON/OD E 2S/JT BLS , STBLS (90)  , DSTA R S (90)  ,0I1DS  (90) ,DBLTS(90)  , 
$OT?DS  ( 90) 

COBBON/ODE20/JTBLU,STBLO  (90)  , DST  AR  0 (9  0)  , 01 1 DO  (90)  ,DEIT0(90)  , 
*01200  (99) 

COBBON/PPTNVn’E,  NORRPR,OERR  , ITBAX 
COBBON/SPI, YN/YX,WT,OWT,DDWT,  ISETOP,KBID 
COBBON/TFBP1/XCX  .IWALLV 


/ 


5 


: 
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"OUTV  ALFNCE  (RATEM  1 p)  ,DDELDX)  , (VA  LSH 1( 1) ,PFLT) 
I WA  I.L  V=T  W ALL 


■THE  CHANNEL  WIDTH  IS  OBTAINED  VIA  CALL  TO  SPLINE. 
FOUIVALENCTNG  T N PLTCITLY  SETS  D ELT ,H B, UT,U T EQUAL  TO  VALS. 
WH=  DT  FEUS  FR  HEIGHT,  AND,  IF  NOT  SPECIFIED  IS  SET 
TO  AN  AVERAGE  ASPECT  RATIO  OF  8.  XC  IS  THE  LOCATION 
OF  THE  FICTITIOUS  SOURCE  TO  CORRECT  POR  3-D  EFFECTS. 

XCNX  =X'~-  X , THE  DISTANCE  OF  THE  SOURCE  FROM  PRESENT  X. 

-SET  UP  THE  CHANNEL  WIDTH  AS  A FUNCTION  OP  XX. 


,1TBL  = 0 
TSETU P=0 

IMIWALL.EO.pGO  TO  25 

IMXX.GT.O)  , THEN  THIS  IS  A CONTINUATION  OF  BL  ROUTINE 

FOR  WHICH  H >HSFP. 

IF  (XX. GT. 0.0) JTBL= JTRLS-1 
X = TX 

TNAX=SW(NPC*1) 

KNID=JSTRTS+1 

CALI.  SPLINE  (SW,  WT,  DWI , DDBI ,D S , JSTP TS, JE NDS,  NDID.4) 

GO  TO  3P 

25  I F (XX. GT. 0.0) JTPL=JTBLU- 1 
X=XX 

XNAX  = SWU  (N-NLC) 

CALL  SPLTNE  (S  WU  , WIU  , DW  I,  DDW I , DS  , JST  RTU,  J EN  DU,  N DIN  , 4) 

30  CONTINUE 


INITIALIZE  COUNTERS  AND  COHPUTE  THE  START  VALUES  OF  DELT,UB, UT. 

START  VALUE  FOP  UT  WAS  READ  IN  AND  PASSED  THRU  CONHON/DER2/ 

N L 00  P = 0 
IRUNGR=1 
DX=DE  LST 
I EX  TT  = 0 
UIREF =UT 1 

0=UTRFF* ( W 1-N  BL* DELST 1 ) 

C PCOEP=1 .0/ (1 . 0-N BL*DE LST 1/W 1) **2 
WRIT*'  (8,040)  TJIPFP,  Q,CPCOEF 
«H=AS*W1 
X C NX  = 1 . E4 

IF  THIS  IS  NEW  RUN,  CALL  START. 

WT=UT 

IP(XX.GT.O.O)  GO  TO  40 
nWT=n. o 
CALL  START 
40  DO  50  J = 1 , 4 

VALS(J)=VALSH1 (J) 

IF(XX.F0.0.0)PATEM1(J)  =0,0 
50  RATES  (.1)  = RATED  1 (J) 

INITIALIZE  LAG  PARAHS  WTTH  EQHTLIBRIUD  VALUES. 

TF  (XX . GT. 0,0) GO  TO  60 
CALL  TA  UN  A X (DUNN  V) 

XO=X 

TAUNO=TAUHR0 
T AUD=T AUH  FQ 

SP  WRTTtf  (6, 400) 

PRTNT  INITIAL  VALUES. 

GO  TO  105 


I BEGIN  D AIN  LOOP 

100  NLOOP=NI OOP* 1 


n n 


J 


Z STORE  CURRENT  VALUES. 

DO  110  .1  = 1,4 
VALSHI(J)  =V ALS  (J) 

HO  R ATE  HI  (J)  = RATES  (J) 

X Y=X 

C STORE  THE  LAG  PAPEHBTERS. 

XO=T 

T AOHO=T  AOH 

IP  ((UB-HI)  FfJT.LE.O  ,0)GO  TO  10.1 

C CHANGE  THE  SIGN  OF  OT  TO  REHOVE  THF  DOU BLE-V ALU EDN ESS  OF  UT . 

UT=-nT 
VALS  (1)  =UT 
WT=T1T 

CALL  BLVALU 
101  CONTINUE 

C PRINT  CURRENT  VALUES. 

TP  (HOD (NLOOP,IPR)  . NE.O) GO  TO  15rt 
1 OS  JTPL=JTBL  *1 

HSFP=  1 . ♦ 1 ./  (1  .-  DELS^/DELT) 

DODX=10.O*TAUH/UT**2 
CP=1 . 0-  (UI/UIREF) 

WRITE  (6,010) XX,DELST,H,HSEP,CP, DELT,UB,UT,U I, CFD2, DQCX 
TF(IWALL.  EO.  1)GO  TO  110 
STBLS  (JT  BL) = XX 
DSTARS  (JT  BL)  = DEL  ST 
0 I IDS  (JTBl)  =UI 
DELTS (JTBL) =DELT 
GO  TO  10  0 
110  STBLU  (JTBL)  =XX 

nSTARU  (JTRL)  =DELST 
U 1 1 DU  (JTBL)  =U  I 
DPLTU  (JTBL)  =DELT 
140  CONTINUE 

TE(IEXIT. EQ. 1 ) GO  TO  260 
160  CALL  ADANS  (DX ,4 , DEP SI , IRU NGE) 

IF(X.  LT.XHAX)  GO  TO  100 

EX  TT  VALUE  CALCULATIONS . 

OBTAIN  VALUES  AT  X=XHAX  BY  EXTRAPOLATION. 

I EXIT= 1 
X X = XH AX 
DO  16  0 .1=1,1 
J J*.1T  BL- 1 ♦ J 
VP(J)  =S"'BTS(.1J) 

Y P (J)  = DREARS  (.1.1) 

UP(.I)  =DELTS  (J  J) 

160  7.  P(.l)  =UT1DS  (.1.1) 

D E LST  = YT  NT  (XP, YP, XI  AX) 

U I=Y INT f X P ,2P  , XI  AX) 

DELT  = YTN-  (XP,UP,XHAX) 

TF  LAS'”  TWO  VALUES  OF  y ARE  VERY  CLOSE  TOGETHER,  THEN  CAN  HAVE 

C PPORLFNS  WITH  S PLT NE-ELI H IN  AT E LAST  BUT  ONE  POINT. 

T”(XN  AX- STBLS  (JTBL-1)  . LT . 0 . 0 0 6*  ( ST  BI.S  (JTBL-1)  -STBLS  (JTBL- 2)  ) ) 

$ JTBL= JTRL-1 
JTPLS=.7TBL  ♦ 1 
GO  TO  IDS 
260  CONTINUE 

900  *ORN  A T ( 1 H 0 , • X DSTAH  H HSEP  CP  DE 

*LT  UB  UT  UT  C P/2  DODX/UI '/) 

010  FO°NATnP1').S,UX,  ?F7.1,2X,F6.  1, 4F 1 2 . 5 ,P1 2 . 6 , FI  2 . S) 

qu*  F0P<|(1»(1«  ,•  REFrpENcr  OU  A NT  TT  T ES,  V EL  OCITY  = • , F 1 0 . 6,  • VOLUfl  E FLOWRAT 
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**=• ,pip.s , «cp  nultiplier=  • , fit . s> 

P ETUP  N 
PND 

S UPPOUTI  N E TPLSTL(IWALL) 

--CALCULATE  DST  ARS  ,01  IDS  FOR  A TBL  AND  1 -D  COPE  IN  SIMULTANEOUS 

ttepation,  assuming  linear  vel  profile  across  channel. 

TWALL=P  IS  LOWER  WALL.  TW  A LI.=  1 IS  UPPER. 

EXTERNAL  DERSTL 

REAL  K AP,VALSM  1 (4)  ,RATEM1  f4) 

REAL  XP(3)  ,YP(3)  ,ZP(3)  ,UP<3) 

INTEGER  KNID,JSTAPT,JEND,JTBt,IRUNGE,IEXIT 
COMMON/A DAM1/X  ,VALS (4) , RATPS  (4)  , PATE  (4, 8) 

COHNON/DEP 1/DDELDX ,DURDX ,DUTDX,DUI DX 

C0MM0N/DER2/DFLT, UB , UT , U I, VT , VB , UDUI , TA UM,H , THETA , DEIST ,CFD2 , 
*VTSCOS,NPL 

COMMON/GEOM 1/ XD , W 1 ,TH, WI DT H, X 1 , B 1 , S I NTH ,COSTH ,TWOTHR  ,TWOTHl , 
«SINTH 1 , COST HI  , AS , W H ,XC  1 , X 2,X  N AX , XDE 
COMMON/OEOM2/N, NR, NL, NU , N M 1 , N LC, NR C 
'-OMNON/TN  TT  AL/DELST  1 , Hi,  U II,  I PR  1 
COMMON/LAG/XO ,TAUMO ,TAUME0 

CONNON/f.  I N E AR/W  DI  F (90)  ,DU2D(90)  ,DDU2D  (90)  , W ND  , D WED  X , DEEP  , DUE  CX 
COMMON/ODF1  S/.1STRTS,  JENPS,NDIM,  SW(90)  ,W  1(90)  , DWI  (90)  ,DDWI  (90)  , 
*DS  (40) 

CONNON/ODE1U/JSTRTU, JENDU.SWU  (90)  ,WIU  (90) 

COM"ON/ODE?S/JTBLS, STBLS (9C)  , DST APS  (9  0)  , U 1 1 DS  (90)  , DELTS  (90)  , 

«U T2DS  (90) 

CONNON/ODE2  U/.1TPLU  , STB  LU  (90)  , DSTAPU  (9  0)  , U1 1 DO  (90 ) , DELTU  (90)  , 
SUT20U (90) 

COMMON/PR TNT/I PP, NORN  PR, CPEPOR, TTN  AX 
C ONE ON/S PLYN/XX,WT,DWT,DDWT, I SETUP, KM  ID 
COHNON/TE  NP1/XC.IWALLV 

EQUTV  A PENCE  (PATEN1  (1)  , DDELDX)  , (VALSN1(1)  , DELT) 

IWALLV=IW ALL 

--THE  CHANNEL  WIDTH  IS  03TAINED  VIA  CALL  TO  SPLINE. 

TOUT  VALE  NCI NG  IMPLICITLY  SETS  D ELT ,U B, UT,U I EQUAL  TO  VALS. 

W H = P T F P 0 S FR  HEIGHT,  AND,  TP  NOT  SPECIPIED  IS  S FT 
TO  AN  AVERAGE  ASPFCT  RATIO  OF  B.  XC  IS  THF  LOCATION 
OP  THF  FICTITIOUS  SOURCE  TO  CORRECT  FOR  3-D  EPPECTS. 

XC  IS  THF  DISTANCE  OP  THE  SOURCE  FROM  THE  ORIGIN. 

W TH  CONTAINS  UT2DU,  WT  CONTAINS  (WDIP-DSTARU) 

,TTBL  = 0 
TSETHP=0 

--TF(XT.GT.O) , THEN  THIS  IS  A CONTINUATION  OF  BL  ROUTINE 
FOP  WHICH  CP>0 . 3 . 

TE  (XX. GT,  0.0)  ,JTBL=  JTBLS-1 
X=XX 

TNAT=SW (N  RC* 1 ) 

KMTO  = .YSTRTS»1 

CALL  SPLINF(SW,WI,DWI,DDWI,DS, JSTPTS, JENDS ,NDIN , 4) 

TSFTUP=0 

call  spline  (SW,  Wit,  DU2D,DDU2D,nS,JSTRTS,  JE  NDS  , N DI  N , 4) 

--INITIALIZE  COUNTERS  AND  COMPUTE  THE  START  VALUES  OF  DELT,UB , UT. 
STA  PT  VALUE  FOP  ni  WAS  READ  IN  AND  PASSED  THRU  COMHON/DER  2/ 
NT.OOP=C 
T RHNGE=1 
DX  = DEI.ST 
TEXIT*0 
U IREE=UI 1 
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0=UTRRF* (W1-N  BL*DELST 1) 

C PC0EF=1 . 0/(1 . 0-HBL*DBLSTl/W  1)  **  2 
WRITE  (6,940) UIREF ,Q .CPCOEF 
WH=AS*W1 

C TP  THIS  IS  HEW  RON,  CALL  START. 

WT*UI 

IP(XX.GT.O.O)  GO  TO  40 

owt=o.o 

CALL  START 
40  00  50  J=1 , 4 

V ALS(J)=VALSH1  (J) 

IF  (XX.  EQ.0.0)  RATEH1  (J)  =0.0 
50  RATES  (J)  =RATEH1  (J) 

C INITIALIZE  LAG  PARAHS  WITH  EQUILIBRIUH  VALOFS. 

IF (XX.GT. 0. 0) GO  TO  60 
CALL  TAOHAX  (DUHH Y) 

XO=X 

TAOHO=TAOHEQ 
TAUH=TAUH  FQ 
60  WRITE  (6,900) 

CFD2=  (OT/OI)  **2 

C PRINT  INITIAL  VALUES. 

GO  TO  105 

BEGIN  HA  I N LOOP  **************** 

100  NLOOP=NLOOP+ 1 

STORE  CURRENT  VALUES. 

DO  110  J= 1 , 4 
VALSH  1 (J)  =V  ALS(J) 

110  RATBH1  (.7)  =RATES(J) 

xx=x 

STORE  THE  LAG  PARENETSRS. 

XO=X 

T AUHO=TAUH 

IP  ( (UB-UI) *UT. LE.0.0) GO  TO  103 

CHANGE  THE  SIGN  OF  UT  TO  REHOVE  THE  DOU BL E- V AL U EDN ESS  OP  IJT. 
U T=-UT 
VALS  (3)  =UT 
WT  = UI 

CALL  BLVALU 
100  CONTINUE 

PRINT  CURRENT  VALUES. 

IE  (HOD  (NLOOP, I PR)  .NE.O) GO  TO  150 
105  .7TBL=JTB LE 1 

HS  EP=  1 . ♦ 1 . / ( 1 . - DEL  ST/D  ELT) 

DQDX=10.0*TAUH/ni**2 
cp=1. 0-  (UI/UTREF)  **2 

WRITE  (6,910)  XX,DELST,H,HSEP,CP,  DELT.UB,  IIT,UI,  CFC2,DQDX 
STBLS  (.1TBL)  =XX 
0 STARS  (JTBL)  = DFLST 
UT1DS  (JT8L)  =0 1 
DELTS  (JTBL) =DELT 
IF(IEXTE. EO. 1) GO  TO  260 
150  CALL  ADAHS  (OX,4,DEPSIL,IRUNGE) 

I E (X.  LT. XH AX) GO  TO  100 

EX TT  VALUE  CALCULATIONS. 

OBTAIN  VALUES  AT  X=XHAX  BY  EXTRAPOLATION. 

TEX  TT  = 1 
X X=X  H AY 

DO  160  j=1,7 
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J J=JT  RL-  1 +J 
XP  ( I)  =ST R LS  (JJ) 
vn  (j)  =BSTARS  (JJ) 

H[»( J)  = DFT.  TS  (JJ) 

If'11  7.  P (J)  =1111  DS  (JJ) 

OELST=YTNT  (XP,YP,XMAX) 

UT=YTNT(YP,7P,XMAX) 

DELT=YTNT  (XP  , UP  , X 1 A X ) 

TF  LAST  TWO  VALUES  OF  X ARE  VEPY  CLOSE  TOGETHER,  THEN  CAN  HAVE 

PROBLEM  S WITH  SPL IN E- EL IM TN ATF  LAST  BUT  ONE  POINT. 
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